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We generalize and extend the stochastic path integral formalism and action principle for con¬ 
tinuous quantum measurement introduced in [A. Chantasri, J. Dressel and A. N. Jordan, Phys. 
Rev. A 88, 042110 (2013)], where the optimal dynamics, such as the most likely paths, is obtained 
by extremizing the action of the stochastic path integral. In this work, we apply exact functional 
methods as well as develop a perturbative approach to investigate the statistical behaviour of con¬ 
tinuous quantum measurement. Examples are given for the qubit case. For qubit measurement with 
zero qubit Hamiltonian, we hnd analytic solutions for average trajectories and their variances while 
conditioning on fixed initial and final states. For qubit measurement with unitary evolution, we 
use the perturbation method to compute expectation values, variances, and multi-time correlation 
functions of qubit trajectories in the short-time regime. Moreover, we consider continuous qubit 
measurement with feedback control, using the action principle to investigate the global dynamics 
of its most likely paths, and Ending that in an ideal case, qubit state stabilization at any desired 
pure state is possible with linear feedback. We also illustrate the power of the functional method 
by computing correlation functions for the qubit trajectories with a feedback loop to stabilize the 
qubit Rabi frequency. 


I. INTRODUCTION 

Continuous measurement of quantum systems [1, 2], 
the study of quantum systems states under the influence 
of observation prolonged in time, has been a topic of con¬ 
siderable activity in recent years. Particularly, for the 
measurement of an individual microscopic system that is 
weakly coupled to measurement apparatus, the system’s 
state and its conditioned evolution in time, the so-called 
diffusive-type quantum state trajectory [3-6], have been 
intensively explored for applications in quantum infor¬ 
mation and quantum control. Some of the active fields 
are quantum state estimation (i.e., estimating the pre¬ 
measurement state of an ensemble of identically prepared 
systems [7]), conditional reversal of measurement [8], and 
the preparation of entangled states [9]. The ability to 
continuously measure quantum systems also opens the 
possibility of feedback control [10], which has also been 
investigated for topics such as the stabilization of coher¬ 
ent oscillations [11-14] and rapid state purification [15]. 

This growing interest in the quantum systems under 
weak continuous measurement has motivated a thorough 
analysis of quantum trajectory statistics. Of notable im¬ 
portance are advances in experiments, such as the mea¬ 
surement of superconducting qubits [16, 17], which has 
allowed the tracking of the trajectories of the quantum 
state with high fidelity in a single measurement run, al¬ 
lowing the statistics of selected subensembles of trajec¬ 
tories to be explored. The authors recently developed 
an action principle [18] over a doubled quantum state 
space, based on a path integral representation of proba¬ 
bility distributions of quantum trajectories. The action 
principle, implemented by extremizing of the stochastic 
path integral’s action, was used to investigate the optimal 
behaviour of the trajectories with arbitrary constraints. 


such as fixing the final boundary condition [19, 20]. The 
stochastic path integral and the optimum likelihood ap¬ 
proach provide a convenient way to investigate statistical 
distributions and globally optimal dynamics of quantum 
state evolution, in addition to the stochastic master equa¬ 
tions describing the quantum trajectories and the Lind- 
blad master equations describing their average evolution 
[3, 21]. 

In this paper, we continue the development of the 
stochastic path integral formalism [18], to further explore 
advantages of having the full joint probability distribu¬ 
tion of quantum trajectories. This includes computing 
statistical averages or expectation values with the abil¬ 
ity to condition on definite quantum states at particular 
times. We present several examples of the formalism in¬ 
cluding a qubit system under the influence of measure¬ 
ment alone, measurement with concurrent unitary dy¬ 
namics, and qubit measurement with feedback control. 
In these examples, the statistics of qubit trajectories are 
computed using developed techniques for the path inte¬ 
gral, such as multi-dimensional Gaussian integrals and 
diagrammatic expansion theory. Moreover, in an exam¬ 
ple of qubit measurement with linear feedback, we utilize 
a phase portrait analysis to investigate the most likely be¬ 
haviour of the system dynamics, revealing a simple and 
practical approach for qubit state stabilization. 

There have been past works on continuous quantum 
measurements with path integrals, so we wish to dis¬ 
cuss how our approach bears both similarities and dif¬ 
ferences to them. An early approach suggested by Feyn¬ 
man [22] and later independently developed by Mensky 
[23] is a restricted path integral: a modified version of 
the Feynman path integral to only sum over paths that 
contribute to a measurement record. Caves, and also 
Barchielli [24-26], constructed similar path integrals by 
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adding coarse-graining (resolution) functions describing 
the effect of the measurement. In these path integral 
approaches, one considers the distribution of the mea¬ 
surement records (it can be derived from the probability 
amplitude, see Appendix A), whereas in our approach, 
we formulate a path integral in quantum state space to 
represent a joint probability distribution of the measure¬ 
ment readouts as well as quantum state trajectories. Wei 
and Nazarov discussed a different approach to continuous 
measurements using the Keldysh path integral technique 
[27]. In Breuer and Petruccione’s path integral [2, 28], 
the notion of the sum over pure state paths in Hilbert 
space is applied, resulting in a different type of doubled 
state space that does not yield an action functional. 

The stochastic path integral formalism and the anal¬ 
ysis of its action are also applied in classical stochastic 
processes. For instance, the formalism is used in studying 
the dynamics and distribution of transmitted electronic 
charge [29], the neural network [30], and the large de¬ 
viations from typical behaviours (rare events) [31]. Our 
approach is similar to the Martin-Siggia-Rose formalism 
[32], which involves adding conjugate fields through the 
Fourier integral form of delta functionals enforcing diffn- 
sive dynamics. Notably, one can think of the quantum 
trajectories on a finite dimensional state space as analo¬ 
gous to classical random processes in configuration space 
of that dimension. 

This paper is organized as follows. In section II, 
we review the stochastic path integral formalism and 
its extremized action equations, for a general finite¬ 
dimensional system with a Markovian setup for weak 
continuous measurement. In section III, a specific mea¬ 
surement setup for a qubit is presented, which is used 
throughout this paper. In section IV, we show that, in 
the case without qubit Hamiltonian, we can perform the 
full path integration directly to get the multi-time corre¬ 
lation functions for the preselected and postselected qubit 
state. In section V, the diagrammatic expansion theory 
is presented as an alternative approximation method for 
computing multi-time correlation functions, with exam¬ 
ples for qubit measurement with Rabi oscillation. In sec¬ 
tion VI, qubit measurement with feedback control and its 
optimal dynamics are investigated using the path integral 
and the action principle approaches. The conclusion can 
be found in section VH. A series of supplementary discus¬ 
sions and some detailed calculations that are not included 
in the main text are presented in the Appendixes. 


II. STOCHASTIC PATH INTEGRAL AND ITS 
OPTIMAL PATHS 

We consider the distribution of quantum state trajec¬ 
tories for continuous quantum measurement. A quan¬ 
tum state trajectory, or simply a quantum trajectory, 
is an evolution of a quantum state in time, conditioned 
on a detector readout realization. This conditional state 
trajectory is also known as a solution of stochastic mas¬ 


ter equations, unravelling master equations in Lindblad 
form. Let us discretize the measurement readout into n 
time points and denote as a measurement real¬ 

ization. Each Tfc is a readout obtained between time tk 
and tk+i = tk + St, and is assumed dependent only on a 
quantum state right before its measurement (Markov as¬ 
sumption). We define a series of quantum states {qk}k=o^ 
written as a d-dimensional parametrized vector q, where 
the components are the expansion coefficients of the den¬ 
sity operator p written in some orthogonal operator basis, 
such as the iV^ -1 generalized Gell-Mann matrices dj of a 
iV-state system [33]. For a two-state system, the matrices 
Oj for j = x,y,z are the Pauli matrices, and q = {x, y, z} 
is a vector in Bloch sphere coordinates. The quantum 
state trajectory can be computed with an update equa¬ 
tion of the form, qk+i = S[qk,rk], taking into account the 
measurement back-action from the measurement readout 
Tfc, and also considering the unitary evolution from the 
measured system’s Hamiltonian. 

Since the distribntion of the measurement readout 
only depends on the quantum state right before the 
measurement in this Markovian approach, we can then 
write the joint probability density function (PDF) of 
all measurement outcomes and state trajectories V( = 
-P({9fc}i) 9 o 5 C) given an initial state qo and a 

set of other constraints C as, 

n-1 

'Pc = B(YlPi^k+i\qk,rk)P{rk\qk)- ( 1 ) 

k=0 

This is a time step product of P{rk\qk), the conditional 
probability distribution for the measurement outcome 
Tk given the system state before the measurement qk, 
and P{qk+i\qk,rk) = S'^iqk+i - €[qk,rk]), the (determin¬ 
istic) conditional probability distribution for the quan¬ 
tum state after the measurement, given the state at the 
previous time step and the measurement readout. The 
prefactor = B(^[{qk},{rk}] in Eq. (1) is a function of 
quantum states and readouts at any times, accounting 
for constraints used in selecting a sub-ensemble of the 
quantum trajectories, such as initial-state and final-state 
conditions. 

The benefit of having the joint PDF Eq. (1) is that it 
contains all the statistical information about the system’s 
evolution under measurement, and it allows us to selec¬ 
tively work with sub-ensembles of quantum trajectories 
simply by adding constraints (or conditions) to the joint 
distribution. Statistical moments can be computed from 
this joint PDF by integrating over its variables. For ex¬ 
ample, an expected value of an arbitrary functional A = 
-4[{gfc},{rfc}] is given by {A)^ = Jd[qk]id[rk]o~^'PcA, 
where we define a notation for the multiple integral, 
/d[qfe]i = /d< 7 i---dq„. Direct integration of these quan¬ 
tities using the joint PDF as in Eq. (1), however, can 
be a challenging task even for a simple qubit measure¬ 
ment problem. As such, we are motivated to write the 
joint PDF in a path integral form with an action (or 
exponent), so we can perform the integration using tech- 
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niques developed in quantum theory such as a diagram¬ 
matic perturbation theory. 

The path integral representation of the joint PDF in 
Eq. (1) can be attained by writing the delta functions for 
the state update 5'^{qk+i - £[qk,rk\) for fc = 0 to n- 1 in 
the Fourier integral form, i.e., 5{q) = (l/27rt) e'^^dp 

for each k and each component of the vector q, and then 
rewrite other terms in exponential forms. The conjugate 
variables for those delta functions are denoted by Pk for 
k = 0 to n - 1. We refer to Ref. [18] and its Appendixes 
for a thorough discussion about this transformation and 
the construction of the path integral. As a result, the 
joint PDF is then given in a path integral form, 

'Pi:=J^Jd[pk]o~^exp{S), (2) 

where the integrals are over all possible paths of {pk}o~^ 
and the action S is defined as, 

S = + l-Pk-{qk+i-£[qk,rk]) + lnP(rfc|qfc)|- 

fc=o I J 

( 3 ) 


We note that is an additional term determined by the 
formation of in Eq. (1), and A/” is a prefactor absorbing 
normalization constants. 

For an example, we consider a sub-ensemble of trajec¬ 
tories that obey conditions on the initial and final states. 
The theoretical analysis of this kind of constraint is pre¬ 
sented in Ref. [18]. The initial state is fixed at go = Qi 
and the final state is at This leads to the con¬ 

straint term = Bg^ g^ = 6‘^{qo - g/)i5‘^(g„ - qp) in 
Eq. (1) and the preselected and postselected joint PDF, 
'Pqi,qF = P{{qk}o,{i"k}o~^,qF\qi), which is Written in a 
path integral form as, 

Pqi,qF = ■^J'ki[pk]_iexp{Sqi^qp), (4) 

where the path integral’s action is. 


Sqj,qi, = - P- 1 - {qo - qi) - Pn - {q-a - qp) 

n-1 f 

+ Y, \ - Pk- (qk+i- £[qk,rk]) + lnP{rk\qk) 
k=0 I 


■ ( 5 ) 


We note that, in Eq. (4) and (5), two additional conju¬ 
gate variables, p i and Pn, are introduced, because we 
have written the delta functions in Bg^ g ^, one for the ini¬ 
tial state and another for the hnal state, in the Fourier 
integral form. 

We can investigate the path integral’s largest contri¬ 
bution by solving for its action’s extrema. Taking the 
variation of the action Eq. (5) over all the variables and 
setting it to zero, we obtain a set of difference equations, 


-gfc+i +£[gfc,rfc] = 0, 


-Pk-i + 


dqk 


{pk - £[qk,rk] + InP(rfelgfc)} = 0, 


d 

— {pk-S[qk,rk] + InP(rfelgfc)} = 0. 


(6a) 

(6b) 

(6c) 


The Hrst, second and third equations are from taking 
derivatives over the conjugate variables pk from fc = 0 to 
n-1, over the state variables g^ from fc = 1 to n - 1, 
and over the measurement readout variables rk from fc = 
0 to n-1, respectively. The derivative over the Hnal 
state variable g„ gives a trivial equation -p„ - Pn-i = 0, 
whereas the derivatives of the action over p i and 
force the boundary conditions go = qi and g„ = qp on 
the solutions of Eqs. (6). We note that, in the case of 
no final boundary condition on the quantum state (i.e., 
= 0 or = Bgj = 5^{qQ - g/)), the derivative of the 
action over g„ instead gives Pn-i = 0, a final value of the 
conjugate variable. 

Interestingly, this extremization of the action in 
Eqs. (6) reproduces the Lagrange multiplier method, an 
optimization strategy that accommodates constraints. 
In our case, the optimized function is the last term 
of Eq. (5), Efe=o InP(rfclgfc), subject to the constraints 
qk+i = £[gfc,rfc] for fc = 0,...,n- I (enforcing the de¬ 
terministic state update) with the Lagrange multipli¬ 
ers po, ...,Pn-i, and the boundary constraints go = qi 
and qn = qp with the Lagrange multipliers p_i and p„. 
Therefore, a solution of the difference equations Eqs. (6) 
is a path that optimizes the log-likelihood of a quan¬ 
tum trajectory, InP(rfcjgfc), subject to the indicated 
constraints. The optimal path can be a local maximum, 
a local minimum, or a saddle point in the constrained 
probability space. For the optimal path that represents 
the local maximum, we call it the most likely path or the 
most probable path. 

The optimal path, a solution of the difference equations 
Eq. (6) and their boundary conditions, can be approxi¬ 
mated by taking a time-continuous limit 5t ->• 0, changing 
the difference equations to a set of ordinary differential 
equations, which can then be solved analytically or nu¬ 
merically. This approximation is applicable because its 
solutions, the optimal readouts and the optimal quantum 
paths, are smooth functions of time. In the case of the 
optimal paths that maximize the log-likelihood of pres¬ 
elected and postselected quantum trajectories (the most 
likely paths between two quantum states), these solutions 
have been experimentally verified with a superconducting 
transmon qubit [19]. 

So far, we have presented the joint PDF and the path 
integrals in the time-discrete form. For a more compact 
representation, we introduce a time-continuous version of 
the discrete stochastic path integral, assuming that the 
system is not intrinsically discrete and a well-defined dif¬ 
fusive limit exists. The joint PDF, for example in Eq. (2), 
is written as, 

Vt;=Af Jd[pk] exp(5) Afjvp exp(5), (7) 

where we have defined a notation for functional integrals, 
e.g., fVp = lim^t^o fd[pk]- The action in Eq. (3) is given 
by, 

5 = ^ dt{- p-{q- C[q, r]) + P[q, r]}, (8) 
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where we have introduced qdt = C[q,r]dt as the time- 
continuous version of the state update equation qk+i = 
and defined JF[q,r] as the linear-order ex¬ 
pansion in time of the log-probability, i.e., P{rk\qk) 
exp{StP[qk,rk] + 0{St^)}. A proportionality factor of 
the latter is absorbed into the normalization factor Af. 
The time dependence of the variables in Eq. (7) and ( 8 ) 
is suppressed for simplicity (e.g., q = q{t)). 

In the continuous version of the action, the state 
update equation q = can be derived from the 

first order expansion in 5t of its discrete form q^+i = 
£-[_qk,'i’k\ using the exact readout probability distribu¬ 
tion P(rfc|< 7 fc). As an alternative, one can consider using 
a stochastic master equation of the quantum state, where 
an ideal white noise ^ is introduced with its single Gaus¬ 
sian probability distribution. The latter substitution is 
valid in an idealized noise limit, which we discuss in more 
detail in section V B 1 (where ltd stochastic equations are 
chosen in the diagrammatic perturbation approach) and 
in Appendix F. The choice of the stochastic state equa¬ 
tions needs to be consistent with the readout (or noise) 
probability distribution used in obtaining the functional 
JF[q,r] in Eq. ( 8 ). 

A careful analysis is needed when considering the ex- 
tremization of the action as described in Eqs. ( 6 ), which, 
in the continuum limit, changes to a set of ordinary differ¬ 
ential equations. We note that one can derive this same 
set of differential equations directly from extremizing the 
continuous-version action in Eq. ( 8 ), if the state update 
equation q = C[q,r] is obtained from the first order ex¬ 
pansion in 6t of its discrete form (this approach to the 
derivation is presented in Ref. [18]). For the action con¬ 
structed using the Ito stochastic equations, its extremiza- 
tion does not give a correct set of differential equations 
describing the optimal paths. However, a stochastic path 
integral formulated using the Ito equations is more ad¬ 
vantageous when computing functional integrals, which 
is presented in section V. 


III. THEORETICAL MODEL: QUBIT 
MEASUREMENT 


So far we have presented the formalism in the gen¬ 
eral context, the quantum measurement with measure¬ 
ment readouts {rk} and corresponding quantum states 
{qk}- In order to show examples and demonstrate how 
to compute interesting quantities using the path integral 
formalism, we choose a particular theoretical model, a 
qubit continuously measured by an apparatus assuming 
a weakly responding (or coupled) detector. This theo¬ 
retical setup has long been a subject of interest in theo¬ 
ries and experiments. Some of the physical realizations 
of this system available with current technology are as 
follows: a single electron in a double quantum dot ca- 
pacitively coupled to a quantum point contact detector 
[34-39], a superconducting qubit dispersively coupled to 
a microwave waveguide cavity [16, 17, 40], and quantum 


optics experiments such as a two-level atom monitored 
with homodyne optical measurement [41]. 

In most of these experimental setups, one can charac¬ 
terize the qubit evolution as governed by the measure¬ 
ment back-action, the qubit unitary evolution, and ex¬ 
tra dephasing due to loss of information or added noise. 
We write a time-discrete change for a qubit density ma¬ 
trix p as pk+i = O-yUst ■M.Tk'lPk^, where we define a mea¬ 
surement operation Mr^[p^ = Mr,,pMl^lPr{Mrt,pM}.i^), 

a unitary operation Ust[p^ = pe^^^*, and an extra 

dephasing operation Oy[p\ (we set h=\). 

The measurement back-action on a qubit state, con¬ 
sidering the diffusive measurement as in Refs. [18, 34, 
35], is described by a measurement operator Mr^ = 
(5f/27rrm)^^"‘exp[-(It(rfe -dz)^/4rm], where Tm is a char¬ 
acteristic measurement time taken to separate the two 
Gaussian distributions by two standard deviations. The 
measurement readout distribution is computed from the 
trace of the measurement operator and the qubit state. 


P{rk\pk) = PT{MrkPkMlj, 
St 


( dt Y 

\ 2TrTm ) 

A—)' 

\ 27rrm / 


r\z ■ 


{1 + Zk)l2 




(l-Zfc)/2, (9) 


with the Bloch sphere coordinates of the qubit state given 

by {xfe,yfc,2fc}- 

We denote a qubit Hamiltonian by H = {el2)dz + 
(-A/ 2 )i 7 a; characterizing the unitary evolution during 
the measurement, and we define the dephasing opera¬ 
tion O^Yp] accounting for additional dephasing on the 
system, such as detection inefficiency and dephasing due 
to the environment. We model the dephasing operation 
as an extra dephasing rate 7 on the off-diagonal elements 
of the qubit density matrix p. By expanding the state up¬ 
date equation = O^Ust AAr^'lPk] to first order in St, 
and taking a continuum limit St 0, we obtain a set of 
differential equations, 


X =-^X - ey - X zrlTm, (10a) 

y = -^y + ex + /\z-yzr/Tjn, (10b) 

z =-Ay + {1-z'^)rlTm, ( 10 c) 


which turns out to be analogous to the stochastic master 
equation in Stratonovich interpretation as mentioned in 
Ref. [34, 35]. 

Using the state update equations in Eqs. (10) and the 
readout distribution in Eq. (9), we have the qubit action 
in continuous form (as in Eq. ( 8 )), 

rT 

Sqb = B^ + dt{-px{x + "fx + ey + xzr/Tm,) 

-Py{y + ^y-ex- Az + yzrlxm) 
-Pz{z+ Ay-{l-z‘^) r/Tm) 

-{r'^ -2rz + l)/2Tm}, (11) 
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where we have approximated the logarithm of the readout 
distribution as \nP{rk\zk) ~ -{5tl2Tm){r\ - 2rkZk + 1) + 
(1/2) \n{St/2TTTjn) + 0{St'^), omitting the second term and 
its higher order expansion because they can be absorbed 
into the prefactor Af. We note that the extremization of 
this action and its most likely paths with fixed initial and 
final states are presented in more detail in Ref. [18]. 

In the remainder of this paper, we will focus on apply¬ 
ing the path integral formalism to this particular qubit 
measurement system, though in different regimes of the 
Hamiltonian parameters. For example, in the next sec¬ 
tion, we study the “plain” qubit measurement where we 
assume that there is no qubit Hamiltonian, i.e., e = A = 0, 
and the qubit’s evolution comes only from the measure¬ 
ment and the dephasing mechanism. 


IV. EXPANSION AROUND THE OPTIMAL 
PATH IN THE PLAIN MEASUREMENT CASE 

In this section, we present examples as to how we can 
use the path integral and its optimal path to compute sta¬ 
tistical moments of the measured qubit trajectories. As 
it turns out, analytic solutions are possible for the case 
of quantum non-demolition measurement (A = 0, where 
the system Hamiltonian commutes with total system- 
detector Hamiltonian) with a constraint on the initial 
and final states of the quantum trajectories. In this case, 
the path integral can be written only in the ^-coordinate, 
i.e., q = z, because the evolution of z is independent of 
the other two coordinates, x and y. The x and y degrees 
of freedom can be found directly from r and z, so they 
will be dropped from the discussion until section V. 

In the following subsections, we consider the plain mea¬ 
surement case when e = A = 0 for simplicity. We show 
how to compute the path integral using the preselected 
and postselected joint PDF, P({zfc}o, Zi?|z/). 

We first integrate this joint PDF over all intermediate 
variables to get the probability density function P{zf\zi) 
of arriving at the final state zp after time T, given the 
initial state z/. Then, we show how to generalize the 
integration procedures to compute statistical quantities 
such as averages, variances, and correlation functions of 
the qubit trajectories in this simplified case. We note 
here that even though we present the derivation in the 
plain measurement case, the result should still be valid 
for the case when e 0. In the latter case, the evolution 
of the X and y coordinates of the qubit will be determinis¬ 
tically oscillating with the frequency e, without affecting 
the evolution of the z coordinate. 


A. Probability density function of a final state 
given an initial state 


Pi{zk}oi{rk}o ^,zf\zi) of the qubit trajectories {zk}o 
and measurement outcomes {rk}o~^ with fixed bound¬ 
ary states, and then integrate over all variables except 
the final state zp, 

P{zf\zi) = fd[zk]^d[rk]r^V.F.,, (12) 

where, as before, we have defined the multi-variable in¬ 
tegral, e.g., /d[zfc]o = /dzo---dz„. 

From the discussion in section H, the joint PDF is given 
by a product over the sliced time variables, 'Pz,,zp = 
6{zo - z/)S(zn - zp) nk=o ^(zk+ilzk,rk) P(rklzk)- The 
term describing the state update equation is a delta func¬ 
tion, 

PizkPi\zk,rk) = ZkSinh^^r 

with its argument derived from the z-coordinate of the 
density matrix equation pk+i = OryUstMr^.[pk], whereas, 
the probability distribution for the readout Eq. (9) is 
expressed in this form, 

P{rk\zk) exp|^(rfc-2rfcZfc + I)|, (14) 

neglecting higher orders in 5t. Substituting these two 
expressions, Eqs. (13) and (14), into the joint PDF, we 
get 

'Pzi,zp = S(Zn - Zf)S(zo - Zl) | H <5(^/c+l"-)| 

xexp]-^-— {rl-2rkZk + l)\, (15) 

I fc=0 J 

where we have used i5(zfc+i-") to represent the delta func¬ 
tion in Eq. (13). 

We then integrate out the measurement readouts by 
transforming the delta function of z into a delta function 
in the readout variable, S(zk+i---) = tanh'^ Zk+i + 

^ tanh~^ Zk)Tml{i - Zk+i)^St, using the transformation 
relation S[x-f{y)] = T,^S[y-y^]/\^yf{y)\y.^f-l(^^), where 
the sum extends over all solutions of yi = f~^{x). Af¬ 
ter integrating over the measurement readout and also 
over the boundary state zq, and z„, the joint PDF is 
transformed to a new joint PDF, a function of only the 
intermediate z-variables, 

- (^)’(n 

(16) 

where the action S is. 


In order to compute the probability distribution of 
the final qubit state given the initial state separated 
by time T, we start with the joint PDF Vz^zp = 


s = -5ty |z^ (^fcpi uk? _ 

^01 2 St 2Tm J 

(17) 
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We note that we have simplified the above equation by 
defining a new variable Uk = tanh”^ Zk to write the joint 
PDF in a compact form. This joint PDF of the variable 
is one of the main results shown in this paper. It is 
also important to point out that since we have already 
integrated over the boundary state variable zq and Zn, 
the delta functions for the boundary states have applied 
to the states Zg = zj and Zn = zp in Eqs. (16) and (17). 

The next step is to perform the integration over the in¬ 
termediate variables Zk for k = l,...,n-l. These integrals 
can be done easily by expanding Zk around the optimal 
solution denoted by Uk- We substitute Uk = Uk + Vk for all 
/c’s in the action, and perform integration by parts with 
the vanishing boundaries, tjq = ijn = 0. This substitution 
exactly simplifies to S[u] = 5[m] - T.k=oiVk+i - VkY 
(see Appendix B for more detail). We then change the 
integral measures to drjk = dzfe/(l - zf.) for k = 1,..., n - 1 
and write the probability distribution in Eq. (16) in terms 
of ? 7 -variables. 




2 exp(5 [u]) 


n-1 


xexp]-^ 'ZiVk+i-Vkf 


k=0 


(18) 


where the optimal path Uk is a solution of the extrem- 
ization of the action, which gives Y.k=oi^k+i - ‘^Uk + 
= 0. The solution is Uk = citk + with the 
two constants of integration c\ = ;]^(tanh"^ z/T’-tanh"^ Zj) 
and C 2 = tanh~^ zj. The optimal solution (the most likely 
path) is given by. 


Uk = ^(tanh ^ 


y T-i — fdnVj ^ 


where we wrote tk = k6t. We note that in the 
time-continuum limit, this extremization is equivalent 
to the vanishing functional derivative of the action, 
6S[u]/Su{t) = TmU = 0, which is similar to the Euler- 
Lagrange equation for the classical trajectory of a free 
particle in the u coordinate transformation. 

The integrals in Eq. (12), with the probability distri¬ 
bution in ry-variable Eq. (18), is now in this form, 


P{zf\zi) = jd[r^k]T^P{{r)k]T\zF\zi). (20) 


Eortunately, they are Gaussian integrals in multiple di¬ 
mensions which can be calculated from the matrix in¬ 
tegral, /dJ 7 exp(-|? 7 ^ ■ M ■ T]) = (27r)^/(DetM)^/^, 
where, in our case, f dr/ = fd[rik]i~^ = /d? 7 i"-d? 7 „_i, and 
the matrix M and the vector r] are given by. 


1 

/ 2 -1 0 

...\ 


< ryi ' 

II 

-12-1 

0-12 


, and ry = 


1 

i, ; ; ; 





noting that DetTVf = n{Tmj5t)^ 


After substituting the optimal solution Uk into 5 [m] 
and performing the integrals over rji, ...,? 7 „_i in Eq. (20), 
we obtain the distribution of the final state fixing the 
initial state and the duration of time = T, 


P{zf\zi) 



exp 



(r^ + 1) H— In 
^ ^ 2 



where we have defined f = ^(tanh"^zp’ - tanh”^ Z/) 
as the optimal measurement readout. This solu¬ 
tion is exactly the same as what we would get from 
the change of variables of the probability distribu¬ 
tion, P{zF\zi)dzF = P{rtot\zi)drtot, where rtot = 
{l/n)n:ork = (Tm/T)(tanh ^ zf - tanh ^ z/) is a time- 
average measurement readout. The derivation of the lat¬ 
ter distribution is presented in the methods section of 
Ref. [19]. 


B. Means and variances of qubit trajectories in the 
plain measurement case 


We have derived the probability distribution for the 
measured qubit given the fixed initial and final states, as 
shown in Eq. (16) and (18). In this subsection, we con¬ 
sider computing the qubit trajectory’s statistical quanti¬ 
ties that can be written in this expectation form. 


zp('A)zj - J' d[zk]i 


P{{zk}i '^,zf\zi) 
P{zf\zi) 


\/DetAf 

(27r)"^ 


f d[rjk]r^Ae--^^"-^-^, 


(23) 


where A is an arbitrary functional of Zj at any time tj for 
j = 1,..., n -1 and the matrix M is the same as defined in 
Eq. (21). Note that we have used the notation zp{---)zi 
for a statistical average conditioned on fixed initial and 
final states, z/ and Zf- Since the probability distribution 
in the ? 7 -variables is Gaussian, it is preferable to write 
z in terms of ij, replacing z with z = tanh(M -i- ry), and 
then perform the integration. For example, a conditional 
average of Zj at time tj is given by, 


zp{zj)z, =zp{taAi{uj +'nj))z, 

= tanhuj +ZF{Vj)zj tanh' Uj 

2 ! 


taxda” Uj + (24) 


expanding to second order in ry, where the primes indicate 
derivatives over the u-variable. In the second line, we 
still use the bracket zf{"')zi for the conditional average 
of the variable rj, even if its boundary values are shifted 
to ryo = ry„ = 0. 

The expectation values in terms of ry can be computed 
using the definition of moments in the second line of 
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FIG. 1. (Color online) The comparison between the analytical 
solutions, Eq. (19),(28),(29), and the numerical simulation, 
for the preselected and postselected average (green dashed 
curve), the most likely path (magenta curve), and the variance 
(the two gray thin curves on both sides of the average are the 
average ± standard deviation) of the measured qubit trajecto¬ 
ries. The preselected state (initial state) and the postselected 
state (final state) are zj = 0 and z_f = cos7r/4, respectively. 
The total time is T = 0.6 Tm. The numerical data, plotted 
as gray dots, are analyzed from 5 x 10® trajectories computed 
with time steps of the size St = 0.006 zero qubit Hamilto¬ 
nian, and the final selection tolerance Zf±0.02. The numerical 
data shows excellent agreement with the theoretical solutions. 
The three fluctuating curves are randomly chosen individual 
trajectories. The inset shows the average (dashed green) and 
the most likely path (magenta) on a Bloch sphere plotted in 
the x-z plane. 


Eq. (23). Because it is a multi-dimensional Gaussian in¬ 
tegral, we find the result from Wick’s theorem, 

ZF{VjlVj2"'Vj2m)zi ~ X! 

all possible 
pairings of 
{h J2,...,j2m} 

(25) 

where the left hand side is an expectation value of even 
numbers of r] at times ti,^ 2 : ■■■;^ 2 mj and the right hand 
side is a product of m elements of the inverse matrix 
M~^, summing over all possible pairings between the rj’s 
on the left side. Any other statistical averages with odd 
numbers of ry will vanish. As an example, pairing four 
ry-variables, ZF{'nj'nkmVm)zn yields the sum, + 

The elements of the inverse matrix are of a simple 
form, = {5tlTm)j{n - k)/n for k > j (the 

full matrix is presented in Appendix C). Knowing these 
matrix elements, we can then compute any expectation 
values of the type shown in Eq. (25), such as a two-time 
correlation, 

= ( 26 ) 


for k > j, and statistical moments of even orders of ry, 

»GT,.(2p-1)!!(E)'(i-|)7 (27) 

where we have used the discrete time notation tj = j6t 
and T = tn = nSt, and the double factorial prefactor {2p- 
1)!! = for a positive integer p. This double factorial 
comes from the number of ways of pairing 2p identical 
variables. 

Substituting these quantities into the expectation 
value Eq. (24), we then obtain the preselected and post- 
selected average of Zj, 

ZF {^j )zi ~ tanh Uj - sech^Uj tanh uj , (28) 

Tm\ T! 

keeping terms up to first order of T/r^, which is equiv¬ 
alent to second order of ry because the conditional av¬ 
erage of the second order of ry scales as T/r^, i.e., 
ZF{VjVj)zi ^ TjTm- This solution Eq. (28) is justified 
in a weak-coupling limit » T), however, the full so¬ 
lution to all orders can be computed and is presented 
in Appendix D. Using the same approximation, we can 
also compute another interesting quantity, the variance, 
keeping terms up to first order oiTlr^, 

zpi^Zj) zj. = Zpizj) Z, ~ Zi) 

(^jsech^Uj, (29) 

where its full solution to all orders in T/r^ is also pre¬ 
sented in Appendix D. Finally, we compute a two-time 
correlation in z-coordinate, keeping up to first order in 
TjTra, 

ZF{zjZk)z, Ri — (l- ^'jsech^ujsech^ufc, (30) 

for tk > tj^ which decreases linearly in for a fixed tj. 

We show in Figure 1, the average ZF{zj)zi, its variance, 
and the optimal path (the most likely path) ttj, for the 
preselected and postselected trajectories, compared with 
data from a numerical simulation using the Monte Carlo 
method (see Appendix E for more detail). We generate 
5 X 10® trajectories with a postselection time T = 0.6 
showing an excellent agreement with the analytical solu¬ 
tions. 


V. DIAGRAMMATIC EXPANSION THEORY 

We have shown in the previous section that the pre¬ 
selected and postselected moments (or expected values) 
for quantum trajectory variables, in the plain measure¬ 
ment case, can be computed by expanding the path in¬ 
tegral around its optimal path. Since the action is Gaus¬ 
sian, however, the path integral approach can be useful 
in other cases as well, such as to compute the moments 









and expectation values of qubit trajectories when there 
is no final state condition (state postselection), or with 
non-zero qubit Hamiltonian. In this section, we present 
an alternative method using a diagrammatic perturba¬ 
tion expansion of the action, similar to how Feynman di¬ 
agrams are used in computing path integrals in quantum 
field theory. 

In the following, we start with a brief introduction 
to the perturbation method (a standard method used 
in quantum field theory) and show how one can com¬ 
pute expectation values of system variables from moment 
generating functionals. We will then go on to the ex¬ 
amples, for the case of qubit trajectories with non-zero 
Hamiltonian and no state postselection. We note that 
even though most of the derivations are shown in time- 
continuous form for simplicity, the most straightforward 
derivations are in time-discretized version, and some of 
these are shown in Appendixes G and H. 


A. Brief introduction to the perturbation 
expansion 

Let us suppose that a quantity of interest is an expec¬ 
tation value given in a path integral form, 

{A) = MjvXVXe^A, (31) 

where the quantity A is a functional of a system variable 
A, the integral’s action 5 is a functional of the system 
variable X and its conjugate variable A, and A/” is a con¬ 
stant normalized factor. Although we are writing A and 
A as one-dimensional variables here, a generalization to 
arbitrary dimensional vectors is straightforward. 

To compute the integral above, we write the action as 
a sum of two parts, S = Sp + Sj, one being terms in the 
action that have bilinear forms, (e.g., ~ AA), which we 
call the free action, 

Sp = -J dtdt’X{t)G~^{t,t')X{t'), (32) 

and call the rest of terms in the action the interaction 
action Sj. We then define a free generating functional 
Zp^J, J], a path integral of the free action adding extra 
source terms with functions J and J, 

Zp[J, J] = M 

= exp jy'dMt'J(t)G(t, t')J(t')|, (33) 

where G{t,t') is an inverse of G~^{t,t') in Eq. (32) sat¬ 
isfying fdt"G~^{t,t”)G{t",t') = S{t - t'), and we as¬ 
sume that the Gaussian integrals in the first line pro¬ 
duce a factor Af~^ that leaves no prefactor in the sec¬ 
ond line. Note that these integrals converge when A is 
purely imaginary. Moreover, if the bilinear terms are in¬ 
stead quadratic terms (e.g., ~ A^), the free action will be 


of the form iS^.quad = -ydtdt'X{t)G-\t,t')X{t'), and 
there will be only one source term fdtX{t)J{t) lead¬ 
ing to a different generating functional Zpquad[J, J] = 
expji fdtdf J{t)G{t,t')J{t')]. 

The generating functional Eq. (33) is then used to com¬ 
pute a, free moment defined as = MjVXVXe^’^■■■. 

The free moment of the variable A(t) (or X{t)) at time 
t is simply a functional derivative of the generating func¬ 
tional over the variable J{t) (or J{t)), taking both vari¬ 
ables J and J to be zero at the end. By looking at the 
second line of Eq. (33), one can see that the simplest non¬ 
vanishing free moment is a two-point correlation func¬ 
tion (A(t)A(t'))F = = GitA'y 

Generalizing this to moments of multiple time points, we 
get 


(A(t,jA(t,J-A(t,,_jA(t,,„))F 

(5 J (5 _~ 

5J{tj2^_i) SJftj^^) J 

~ X/ 1=27,1-1 ’ ^^<=2™ 

all pairings between 
variables X and X 


where the summation is for all possible pairings between 
the variables A’s and their conjugate A’s. The number 
of propagators in the summation on the right is equal to 
the number of pairs presented in the expectation bracket 
on the left. From this, one can see that if there is at least 
one unpaired variable, the moment will vanish. 

Using these definitions of the free moments, one can 
then compute a path integral as in Eq. (31) where S = 
Sp + Si and the arbitrary functional A is a function of 
the physical variable A, 


(A[A])EAf 



=(A[A]e 


e‘5^[w’A]ZF[J, J] 
5 i[X.X]^F, 


J=j=0 


(34) 


where in the second line we wrote the interaction action 
Sj and the arbitrary functional A as operators acting on 
the free generating functional Zp[J, J]. 

The perturbation expansion refers to the series expan¬ 
sion of the term exp{iS/[A, A]} in Eq. (34). In some 
cases, infinite series can be summed over, giving an ex¬ 
act analytic result, although in many others, an approx¬ 
imation is needed in order to truncate the series. An 
approximation can be made when there is a small pa¬ 
rameter appearing in any (or all) of the terms in the 
interaction action Sj, where the expansion is straight¬ 
forward, keeping terms up to any desired order of the 
small parameter. Another type of approximation, which 
we present here in more detail, is similar to a semiclas- 
sical expansion in quantum mechanics, keeping terms up 
to any order of a small parameter v that appears as an 
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inverse in front of the action, for example, 

5= + Si[X,X]Y (35) 

where we explicitly write the the free action in a bilinear 
form. To compute an expectation value of a functional A 
using this action, one needs to expand the exponential of 
the interaction action and then evaluate the free moments 
in terms of the propagators. The order of expansion is 
controlled by the parameter v. There is a factor oi Ijv 
for every single term in the expansion of Si, and a fac¬ 
tor of V for every propagator that emerges. This is the 
basic idea behind the loop expansion in quantum theory, 
where v plays the role of the h in the Feynman path in¬ 
tegral. We will discuss more about the loop expansion in 
section VB 2 . 

This loop expansion based on the small parameter v, 
in our quantum measurement case, can be considered as 
a small noise expansion around a saddle point solution (a 
solution that extremizes the action). Moreover, as in the 
Feynman path integral, under the approximation of the 
small parameter v, a saddle point approximation can also 
be applied to estimate a path integral, using an expan¬ 
sion of the action around the saddle point solution. For 
our stochastic path integral, the saddle point approxima¬ 
tion gives an estimation of total probability density for 
the joint probability distribution that the path integral 
represents. 


B. Examples in continuous quantum measurements 


1. Diagrammatic rules for qubit measurement with Rabi 
oscillation (A t 0) and with no final state condition 

The theoretical setup for the qubit with Rabi oscilla¬ 
tion is analogous to the one used in Ref. [19], where the 
qubit Hamiltonian is H = (-Al2)ax- Here we can calcu¬ 
late important quantities such as average trajectories and 
correlation functions, for the case when there is only an 
initial state fixed and not the final state. For simplicity 
of the diagrammatic expansions, we make a white noise 
approximation and variable transformations to the qubit 
system. 

We make an idealized white noise limit on the qubit 
measurement. The variance of the measurement readout 
distribution P{r\z) is assumed to be very broad, which in 
this case means Tm » St, justifying an approximation of 
the two Gaussian distribution in Eq. (9) to a single Gaus¬ 
sian distribution, P{r\z) ^ (( 5 t/ 27 rT™)i/ 2 g-(r-z) 25 i/ 2 r„ 
(see Appendix F for more detail). The measurement 
readout is then approximated to be its mean plus a noise, 
r = z+^Tjn C; where ^ is the Gaussian noise with variance 
St~^, independent of the qubit state 2 . Because the na¬ 
ture of this noise is highly fluctuating, in the derivation 
of the state update equations, we need to keep an expan¬ 
sion up to a second order in St, replacing r'^St^ ~ T^St 
[42]. This leads to the ltd stochastic differential equations 
[34, 35], 

X = -Tx - X , (37a) 

y = -Ty + Az-yz^lrli^, (37b) 

z = -Ay+{l-z^)^lTU\ (37c) 


Now we can apply the perturbation approach to our 
quantum measurement problem. We use the joint prob¬ 
ability density function introduced in Section H in its 
generalized form Vq = P{{qk},{'''k}\Qo,Oj ^ joint PDF 
of the measurement outcomes and the quantum states 
given an arbitrary set of constraints Statistical aver¬ 
ages or expectation values using the joint PDF are then 
given in this form, 

{A)i:= Jd[qk]d[rk]Pc A, 

^*fs^J\fJ'vqVrVp ex.-p{S)A, (36) 

where, in the second line, we have taken the time- 
continuous limit and written the joint PDF in the path 
integral form as Vq = AffVpexp{S), noting that the 
time-continuous action is given by Eq. ( 8 ). 

In the following examples, we present the perturba¬ 
tive expansion approach in computing the statistical mo¬ 
ments of qubit trajectories. We focus on the case when 
the qubit Hamiltonian does not commute with the mea¬ 
surement operators (A t 0) and without a final state 
constraint. 


where x, y, z are Bloch sphere coordinates for the qubit 
and a dephasing rate F is now a total dephasing rate F = 
7 + l/2rm. The white noise ^ has a Gaussian probability 
distribution P(^) = ((5t/27r)^/^ exp(-^^i5t/2). 

Before substituting the state update above into the 
action of the path integral, we can make changes in the 
variables of the system in order to simplify the later per¬ 
turbation process. This is to avoid infinite series related 
to the linear terms in Eqs. (37). We define a new set of 
variables u, v, w where {u, v, w} = Q ^ ■ {x, y, zj and Q is 
a matrix that diagonalizes the linear terms of Eqs. (37), 


/-F 0 

0 ) 

1 / •^i 

0 

0\ 

0 -F 

+A 

0 

II 

A2 

0 

^ 0 -A 

0 j 

' \0 

0 

Asy 


The eigenvalues are Ai = -F, A 2 = -(F + D)/2 and 
A 3 = -(F - D)/2, where we define D = \/r^ - dA^. The 
diagonalizing matrix Q and the transformation of the 
system variables x, y, z are described in the matrix equa¬ 
tion. 


'x^ 


(1 0 

0 \ 

( u' 

y 


n r 4-Q 

r-n 

V 

— 

2A 

V) 


lo 1 

1 1 

\w) 
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Type 

Labels of vertices 

Full forms 

Diagrams 

Type 1 
(initial) 

PuOf PvO^ PwO 

ui Atpu{t)5{t) 


Type 2 

PuUvG PuUW^, PvVVi, 
PvVW^, PwWV^, PnjWW^ 

af^ Atpu{t)u{t)v{t)^{t) 


Type 3 

PuG Pv^, Pwi 

Kif^Atpu{t)i{t) 



TABLE I. We show the 12 terms in the interaction action Eq. (41b) (interaction vertices) classified into three types. The 
second column shows the shorthand labels of the vertices, simply the variables contained in each of interaction terms, while the 
third column presents examples of their full integration forms. In the fourth column, we illustrate the vertices as connecting 
edges. The direction of the arrows shown on solid edges is from the Pu,v,w variables to the u,v,w variables representing the 
propagators Gu,v,-w The wavy lines are for the connection between noise variables. The number of edges (in-coming, out-going, 
and curly edges) around the vertex (small circle) corresponds to the number of connections required for the vertex. 


which gives the following transformation: x = u, y = v {T+ 
n)/2A + re (r - n)/2A, and z = v + w, and the inverse 
transformation: u = x, v = { 2 yA -Tz + 0z)/(20), and 
w = {- 2 yA + Tz + flz)/{ 2 fl). 

The stochastic differential equations (37) with this new 
set of variables are, 

w = Ai w + aw(n + w) ^-t Ki ( 40 a) 

i) = X2V + av{v + w) ^ + K2^, ( 40 b) 

w = X3W + aw{v + w) ^ + K3^, (40c) 

where we have defined a = -l/xm^ and Hi = 0, K 2 = 

(-r + n)i{ 2 TT^), k 3 = (r + n)i{2TU^n). 

Now we are ready to construct the action of the path 
integral, substituting the update equations (40) and the 
log-likelihood function into the action 

Eq. (8). We write the action in two separated terms, 
S = Sp + Si, where 

rT 

Sf = dt{-pu (ii - Xiu)- Py {v - X2V) 

-Puj{w-X3w)-^‘^/2}, (41a) 

rT 

Si = j dt{B + apuu{v + w) ^ + KiPu^ 

+ apv v{v + w)^+ K2 Pv^ 

+ apu,w{v + w)^ + K 3 Pw^}, (41b) 

are the free and interaction actions, respectively. We 
note that the additional term B is not merely a time- 
continuous form of the first term in Eq. (5) (we only con¬ 
sider the initial condition). This is because the actual 
initial condition term -p i • {Qo ~ Qi) can be removed 
simply by integrating over the initial variable Qq, forc¬ 
ing the condition {a;(0),y(0),z(0)} = {xi,yi,zi} to the 
variables at time to = 0. In place of the initial term, 
when we write the free action Eq. (41a) in this form, 
Sp = - fdtdt'X{t)G~^{t,t')X{t'), there will be leftover 


terms which contribute to the term B, resulting in, 

B = J'dtuiPu{t)S{t) + J dtvipy{t)5{t) 

+ J dtwiPn,{t)6{t). (42) 

This can be easily shown with a discretized version of 
the path integral, which is presented in more detail in 
Appendix G. We note also that, in the rest of the paper, 
the time integral is always from t = 0 to t = T, unless 
stated otherwise. 

The propagators (or Green’s functions) are compnted 
from the inverse Green’s functions, which in this par¬ 
ticular case are G^{t,t') = S{t - f) - Xk), where 
k = 1,2,3 (or u,v, w) for the first three terms in Eq. (41a), 
and = 6{t -1') for the noise term. With an identity 
relation, / dt''G~^{t,t'')G{t'',t') = 5{t-t'), the propa¬ 
gators (Green’s functions, or two-point correlation func¬ 
tions) are given by, 

{u{t)pu{t'))F = Gu{t,t') = e(t-t')e^'^‘“* \ (43a) 

{v{t)py{t'))F = Gy{t, t') = 0(t - , (43b) 

{w{t)py,{t'))F = Gy,{t,t') = \ (43c) 

{mat'))F = Gdt,t') = S{t-t'), (43d) 

where t>t' for the first three lines. It is important to note 
that 0(t) is a left continuous Heaviside step function, 
i.e., 0 ( 0 ) = 0 and limt^o+ 0 (^) = 1 , causing the two- 
point correlation functions for u,v,w to vanish when t < 
t'. This is a result of the ltd interpretation we chose in 
writing the state update equation in Eq. (40). This can 
also be verified with the discretized version of the path 
integral and is presented in Appendix H. 

To compute quantities such as an expectation value 
{A) = {Ae^‘)F where A is an arbitrary function writ¬ 
ten in the form of Eq. (34), one needs to expand the 
exponential in a power series of Si, and looks for 
terms that contribute to its result. Fortunately, from the 
two-point correlation functions Eq. (43), we know that 
system variables u,v,w can only connect with the their 
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conjugate variables at earlier times, e.g., u{ti) can only 
connect with Pu{t 2 ) if ti > t 2 - Knowing this helps predict 
which terms will contribute to the sum of the expansion. 
For example, to show that (1) = {e^‘)F = 1, we look at 
terms in the interaction action Eq. (41b) and see that all 
of them contain exactly one conjugate variable in each. 
Looking at the propagators in Eqs. (43), we know that it 
is impossible for any higher order terms in the expansion 
of to have all conjugate variables matched with sys¬ 
tem variables at their later times. This is because every 
time we want to find terms with system variables at later 
times, there will be at least one more unmatched conju¬ 
gate variable appearing. Therefore, there are no other 
contributions apart from 1 , resulting in {e^^)p = 1 . 

In order to make the calculation of moments (A) = 
{Ae^^)F more systematic, we develop diagrammatic rules 
to ease the process of keeping track of the perturbative 
expansion. Each term in the expansion of is a com¬ 
bination of 12 additive terms shown in the interaction 
action Eq. (41b), with repeated appearance also possi¬ 
ble. The 12 interaction terms {interaction vertices) are 
shown in Table I with their labeling names and examples 
of their full integral forms, where we use the latter to 
characterize the vertices into three types shown as the 
three rows. At the end, all combinations of the interac¬ 
tion vertices (all additive terms in the expansion of e^‘) 
are then multiplied with the function of interest A be¬ 
fore taking the free expectation The functional 

^ is a function of system variables and noise variables 
(i.e., u,v,w,^) at all times. These variables are called 
ending vertices. Let us use the word combination for an 
individual term combining interaction vertices (from ad¬ 
ditive terms in the expansion e ^‘) and the ending vertices 
(from the function A). One can determine which com¬ 
binations are non-vanishing and compute their values by 
following these diagrammatic rules: 

• A non-vanishing combination of vertices should 
contain equal number of each of the system vari¬ 
ables {u,v,w) and their conjugates {pu,Pv,Pw), and 
there should be even number of the noise variables 

iO- 

• The system variable can only be connected to its 
conjugate at earlier time, while the noise variable ^ 
connects to another noise variable at the same time. 
The connections are presented as ‘edges’, or lines 
joining vertices. The number of edges for each ver¬ 
tex is equal to the number of variables it contains 
(graphical representations are shown in Table I, in 
the column of ‘Diagrams’). A non-vanishing combi¬ 
nation can have any number of vertices, but edges 
have to all be connected. 

• A non-vanishing combination is presented by con¬ 
nected diagrams with the ending vertices (variables 
with their time arguments being in the range of 
t 6 (0,r]) on the most left and their time argu¬ 
ments clearly stated. The time ordering decreases 


Uti Pu 

*—^- 1 

(a) ^ (b) 


FIG. 2. Examples of diagrams, where the labels indicating the 
variables contained in each of the vertices (we label only the 
ones that are necessary). The filled dots represent the end¬ 
ing vertices, while the open dots represent interaction (initial, 
t = 0) vertices, (a) An ending vertex connects with an initial 
vertex (pvo)- (b) Two ending vertices connect with two inter¬ 
action vertices {Pu^ and Pv^- (c) Two ending vertices connect 
with six interaction vertices (from left-right, top-bottom or¬ 
der: PvVW^, Pvi, Pvi, PwWW^, Pwo and Pwo)- The unequal 
length of horizontal edges following the ending vertices Vti 
and Vt 2 indicate that ti >t 2 - 


from left to right, with variables at the same time 
lining up in the same vertical line, and the vertices 
with initial conditions {t = 0 ) on the most right. 

• A result of the free expectation value {{•■■)f) for 
each combination is computed by taking into ac¬ 
count the constants (i.e., ui,vi, ki, ...) and 
integrals attached to each of the interaction ver¬ 
tices, the propagators (the Green’s function), and 
the numbers of possible ways to connect diagrams. 
From the point of view of graphical diagrams, edges 
represent propagators, and points common to more 
than one edge correspond to integrals over time. 
Moreover, if there are m repeated vertices in a com¬ 
bination, it needs to be multiplied by a factor ^ 
(resulting from coefficients in the expansion of e ^‘). 
However, these factors usually cancel with the num¬ 
ber of possible ways to connect the repeated ver¬ 
tices. 

In Figure 2, we show examples of non-vanishing dia¬ 
grams showing connections of one and two ending ver¬ 
tices. For each of the diagrams, we can determine its 
contribution by explicitly writing out its full form and 
computing its integrals. In the following subsections, we 
present some of the calculations using these diagrams and 
rules to estimate statistical averages of the qubit trajec¬ 
tories. 

We note that it is also possible to construct different 
diagrammatic rules based on different formations of the 
path integrals. For example, one can perform integrals 
over the readout noise (e.g., the noise ^), transforming 
the action in Eq. (41a) and (41b) to a new action that 
contains only the system variables u, v, w. In this case, 
the diagrams will be completely different from the one 
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we presented above, but still giving exactly the same re¬ 
sults for the statistical quantities related to the system 
variables. Here, we choose to present the diagrammatic 
rules with explicit noise variables because of the ability 
to compute correlation functions between the noise and 
the system states, as shown in section VB4. 


loops are as follows: (a) if = 1, / = 0, H = 1, L = 0, (b) 
E = 2, I = 1, V = 2, L = 0, {c) E = 2, I = 6, V = 6, L = 1, 
giving the zeroth(zeroth), zeroth(first), and first(second) 
orders of loop(i^), respectively. 


3. Example A: Average quantum trajectory 


2. Small noise approximation and loop expansion 


The diagrammatic rules help in determining non¬ 
vanishing terms, but we need a systematic approach to 
keep track of the order of expansion, especially when 
there are infinitely many ways to construct diagrams. We 
consider the small noise approximation around the saddle 
point solution mentioned in section V A. To obtain the 
action of the form in Eq. (35), we introduce a small pa¬ 
rameter denoted by i' to the noise term of Eqs. (41), con¬ 
trolling the noise variance in the qubit dynamics, leading 
to a modified action, 

rT 

S' = dt{-pu (u - Aim) - {v - \2v) 

-Pwiw- X3w)-^‘^/2iy} + Si, (44) 

leaving all terms in the interaction action intact. We then 
rescale the conjugate variables by replacing Pujv^ 

Pv Pv I V and Pw Pw 1 (as well as changing variables 
in the source term, e.g., Ju,v,w ^ Ju,v,wlv), noting that 
the rescaling of the auxiliary variables does not affect the 
qubit dynamics. Since all terms in the action Eq. (44), 
except for the noise term, have exactly one conjugate 
variable (see Eqs. (41)), the resulting action is then in 
the desired form, 


5-1 


dt[-pu {u - Aim) -py{v- X2V) 
-pyj{w- Asm;) - ^^/2} + 5/ 


(45) 


As an example of how to use the diagrammatic rules to 
compute statistical expectation values, we show a deriva¬ 
tion of {z{t)), an average of the z-coordinate of the quan¬ 
tum trajectories. Erom the variable transformation in 
Eq. (39), we know that z{t) = v{t) + w{t), therefore we 
can write the expectation value in terms of the free mo¬ 
ments and then compute the diagrams from the vertices 
in Table I, 

{z{t)) ={v{t)e^^)F + {w{t)e^^)F, 

• < o •—<—o 

^ Vt PvO . Wt PwO 

= Vijdt' Py{t')S{t')^^ 

+ ^w{t) Wijdt' p^{t')S{t')^ 

= dt'Gv{t,t')S{t') + wij' dt'Gw{t,t')S{t'), 

= where t > 0, (46) 

where the Green’s functions are from Eqs. (43). Erom 
the types of vertices in Table I, there is only one possi¬ 
bility connecting the ending vertices v{t) and w{t) with 
PvO and PwO, respectively. This is because choosing other 
vertices will lead to at least one unmatched noise vari¬ 
able which, when trying to pair it with another vertex 
containing will result in more unmatched conjugate 
variables. 

As a result, the expectation value Eq. (46) is explic¬ 
itly found by transforming back to the system variables 
x,y,z, 


where Pu,v,w are the rescaled conjugate variables. 

Eollowing the discussion in section V A and the dia¬ 
grammatic rules presented in the previous section, the 
order of expansion can be determined by counting the 
number of edges and vertices of a diagram. Each vertex 
in a diagram contributes a factor of 1 /m to its result, while 
each edge (which represents a propagator) contributes a 
factor of M. Therefore, each term (diagram) in the ex¬ 
pansion carries a factor of , where E, /, and V 

are numbers of its external edges (edges that connect 
with ending vertices), internal edges, and vertices, re¬ 
spectively. This corresponds to the loop expansion where 
one instead counts the number of loops from L = I-V+1, 
which leads to another way of writing the order of the ex¬ 
pansion , noting that the number of external edges 

E is fixed for each set of ending vertices. 

We show how to count the power of v of the diagrams 
shown in Eigure 2. The numbers of edges, vertices, and 


{z{t)) = e I \^zi cosh—+ - - -sinhyj, (47) 

which is exactly the same as the solution we would get 
from averaging over all possible noise realizations in the 
ltd stochastic master equations in Eqs. (37) and solv¬ 
ing for an average of z{t). The averages of x and y 
can be computed in the similar way. We note that the 
exact solution for the average trajectory involves only 
the diagrams with zeroth order of the small parameter m 
iE + I-V = 0). 

Particularly, one can show that the average trajec¬ 
tory coincides with a saddle point solution of the action 
Eq. (44) when m -> 0. The vanishing functional deriva¬ 
tives of the action over the variables Pu,v,w gives the ex¬ 
act same equations Eqs. (40) (though, the noise variable 
^ is no longer a stochastic function), while the vanishing 
functional derivatives over m , v, w leads to the differential 
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equations for conjugate variables, 


Pu = -Aip„ - a{v + w)pu^, (48a) 

Pv = -^2Pv - a{upn + (2u + w)py + wp^}^, (48b) 

Pw = -Asp™ - a{upu + vpv + (u + 2w)pw}^, (48c) 

where the noise variable is now a solution of extremizing 
the action over ^ giving, 

i2{a{v + w){upu + vpy + wpu,) 

+ KiPu + K2 Pv + KsPyu}. (49) 

As the noise parameter i/ decreases to zero, as does the 
noise term the ordinary differential equations of u, v, w 
are uncoupled from Eqs. (48) and their solution (a saddle 
point solution) is then the same as the average trajectory, 
the solution of Eqs. (40) when ^ = 0. 

4- Example B: Correlation function between system 
variables and noise variables 

Another example is a correlation function between the 
system variable 2 and the noise variable given by 

+ {w{ti)^{t 2 )). Since there 
are now two ending vertices, the situation is more com¬ 
plicated and there are infinitely many ways to connect 
the diagrams. However, we can compute terms to the 
lowest order of the loop expansion, which in this case is 
the zeroth loop (tree-level diagrams). There are in total 
six contributions to the correlation function, as shown in 
the following. 



= K2j ,t2) 

+ (Xv'j J' dt Gy {tl,t ')G^ (t , t2')Gy (t , 0)^ + •", 

= + aviwie^^*^) 

+ + avjwie ^^*^) , (50) 

where the superscript ( 0 ) in the first line indicates the 
number of loops in the expansion, and, in the third line, 
we show only the integral form of the first two dia¬ 
grams. This result involves diagrams with first order of 
the small noise parameter i/. The definitions of the pa¬ 
rameters such as t 6 /,K 2 ,A 2 are defined in the discussion 
of Eqs. (38)-(39). 

We note that this result Eq. (50) is valid for ti > t 2 
and it is zero otherwise, because of the special properties 


of the left continuous Heaviside step function 0(t). This 
can be interpreted as the noise being correlated with the 
qubit state at later times but not with the state at earlier 
times. This property of the system variable and noise the 
correlation function is also mentioned in Ref. [43]. 


5. Example C: Correlation function between two system 
variables 


Let us next consider correlation functions between sys¬ 
tem variables at two points in time, such as (z{ti)z{t 2 )) 
and {y{ti)z{t 2 ))- Both of them can be written in terms 
of correlation functions between the transformed sys¬ 
tem variables v,w, for example, the z-z-correlation is, 
{z{ti)z{t2)) = {v{ti)v{t2)) + {v{ti)w{t2)) + {w{ti)v{t2)) + 
{w{ti)w{t 2 ))■ As in the previous example, we keep terms 
up to the lowest order of loops, which is the zeroth order. 

Each term of the correlation functions, for example 
{v{t\)v{t 2 )), involves 10 different diagrams, thus there 
are in total 4x 10 diagrams for computing the z-z correla¬ 
tion function. We present in the following three samples 
of diagrams and their integral forms for the correlation 
function {v{ti)v{t 2 )), 


(v(fi)u(t2))^°^ = {v{ti)v{t2)e^')^P, 


Vti PyO 

• - i -O 


• < o 

^^2 PvO 


+ 



Vt 2 Pv 


+ 



+ (7 more similar diagrams), 

= (v{ti)){v{t2)) 

+ kIJ dt'dt"Gy{ti,t’)Gy{t2,t")G(,{t',t”) 

+ aK2v'j J dt' dt” {Gy{ti,t')Gy{t2,t'') 

G^{t',t'')Gy{t',d)Gy{t',0)} + -, (51) 


where, apart from the first trivial diagram, each initial 
vertex {vt^ or vt^) can connect to three possible interac¬ 
tion vertices {py^, pyVV^, or pyvw^). As a result, there 
are total of 9 + 1 possible ways, and they are presented 
in Appendix I (all diagrams except the first one are of 
the first order of v). The similar kind of calculation is 
applied to the other correlation functions {v{ti)w{t 2 )), 
(w{ti)v{t 2 )) and (w{ti)w{t 2 )). 

The correlation function between y and z can be calcu¬ 
lated in a similar way, with the same types of diagrams, 
but with extra prefactors /3i and (32, 


{y{ti)z{t2)) = I3i{v{ti)v{t2)) + I3i{v{ti)w{t2)) 

+ I32{w{ti)v{t2)) + P2{w{ti)w{t2)), 

where the coefficients /3i = (T + H)/ 2 A and /?2 = (L - 
H)/ 2 A come from the transformation of variables shown 
in Eq. (39). 
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FIG. 3. (Color online) The variances and covariances for the 2 -coordinate of the qubit trajectories, comparing the theoretical 
solutions computed up to fourth order of vertices and numerical simulations (based on 10® trajectories). The qubit Rabi 
frequency used in these cases is A = 207rT,j,^, and the total time is T = 0.5 Tm- The regime presented in the upper row of panels 
is the case of efficient measurement with no extra decoherence (i.e., 7 = 0 , ry = l/( 2 Tmr) = 1 ), while the regime presented in 
the lower row case is with low efficiency 7 = 0.02. The initial state used here is {xi,yi,zj) = (0,1,0). (a), (d) The simulated 
individual trajectories for both regimes, (b), (e) The theoretical variances {z{t)^)-{z{t))^ (red curves) along with the numerical 
data (gray dots). In the panel (e), the variance increases from zero at the initial point t = 0 and then converges to a value of 
(r 2 +A 2 )/( 2 rAV„,) «(0 .023. (c), (f) The covariances between 2 at time ti = ta = 0.113Tm and time t 2 = ta + t (solid magenta 
curve), and between 2 at time ti = tb = 0.339 and time t2 = tb + T (dashed blue curve), are presented along with results from 
the numerical simulation (gray dots). 


6. Discussion and comparison with numerical simulation 

We compare the theoretical results, specifically the 2-2 
correlation computed in the previous section, with nu¬ 
merically simulated quantum trajectories. Although, the 
solution of (2(1:1)2(12)) is not explicitly presented here in 
this paper, as it is too lengthy, it can be computed in a 
similar way as shown in Eq. (51) and also in Appendix I. 

We show in Figure 3(a) and 3(d) sampled individ¬ 
ual trajectories generated from the Monte Carlo method 
(Appendix E), in the regime where T < Tm- The theo¬ 
retical variances are computed from the 2-2 correlation 
function, and they are shown in panels (b) and (e) along 
with the variances from the numerical trajectories. In 
the top panel, for efficient detection (ry e l/(2rmr) = 1), 
the numerical variance increases in time and gets satu¬ 
rated (not shown) at around a value of 0.5 for a long 
enough time, indicating that the qubit states at the later 
times are distributed throughout the perimeter of the y-z 
plane of the Bloch sphere. The variance from the tree- 
level diagram approximation fails to exactly capture the 
long-time behaviour, as we can see that the discrepancy 
in the panel (b) starts to grow as time increases. How¬ 
ever, for the inefficient detection case, rj = 0.02 in Fig¬ 
ure 3(e), we can see that the theoretical approximation 
can explain the behaviour quite well, predicting the satu¬ 


rated variance at a value of (r^ +A^)/(2rA^rm) « 0.023. 
To compute this quantity theoretically, one takes a limit 
t 00 of the calculated variance {z{t)‘^) - {z{t))‘^. 

For the correlation function at two different times, we 
define a covariance as Cov[z{ti)z{t 2 )^ = {z{ti)z{t 2 )) - 
{z{ti)){z{t 2 )), and show in Figure 3(c) and 3(f) its nu¬ 
merical and theoretical comparisons. We plot the covari¬ 
ance for ti = tafi and ^2 = ta,b + t where r has its value 
ranged from 0 to T-ta,b- For the efficient detection case, 
shown in the panel (c), the agreement between the the¬ 
ory and the simulation is better for the short time case, 
ta = O.llSrm. For the inefficient detection case, panel 
(f), the agreement is excellent for both ta and t^. 

Since our theoretical results are derived with the small 
noise expansion around the saddle point solution (or the 
average solution, see section VB3) using v as an expan¬ 
sion tracking parameter, the diagrammatic approxima¬ 
tion method works well whenever the qubit trajectories 
are narrowly distributed around its average. This hap¬ 
pens in the short-time regime when the diffusion is still 
small from the initial state, as well as in the regime when 
the dynamics is strongly suppressed by the dephasing 
mechanism (when the detection efficiency is low). We 
note that for the long-time limit, one needs to com¬ 
pute higher order terms which contain more complex di¬ 
agrams. However, there are other approaches, such as 
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using the qubit master equation to approximate the cor¬ 
relation functions in the stationary limit and calculate 
spectral densities of the measurement readouts. These 
can be found in the works on continuous measurement of 
mesoscopic electronics such as in Refs. [37, 43-45]. Com¬ 
paring to these approaches, our results give much more 
accurate correlation functions at short times. 


VI. STOCHASTIC PATH INTEGRAL FOR 
CONTINUOUS MEASUREMENT WITH 
FEEDBACK 

We can extend the discussion of the stochastic path 
integral formalism to its application to the system un¬ 
der continuous measurement with feedback control. The 
feedback loop consists of getting information about the 
system state via the measurement, and feeding back a 
control signal to the system in order to alter the state as 
desired. One of the most intuitive models of the feedback 
loop is that the system Hamiltonian changes as a function 
of the system state, which then can be written as a func¬ 
tion of the most recent measurement readout. For exam¬ 
ple, in our qubit measurement case, taking into account 
additional time delay r^, the Hamiltonian at any time t 
can be written as a function of the measurement readout 
in the past, Hfb{rt-r^,t) = €{rt-rd)^zl‘2. - A{rt-Td)^xl‘2., 
where the parameters e and A, as defined in Section HI, 
are now functions of the measurement readout at 

time t-Td- 

We consider an ideal case with instantaneous feedback 
Td ^ 0, i.e., the measurement readout at time t imme¬ 
diately changes the system parameters which are used 
in computing the state update at the time. This way, 
the formulation of the stochastic path integral as time- 
local, its action’s extremization equations, and the dia¬ 
grammatic expansion we presented so far are perfectly 
applicable. The only modification needed is to treat the 
readout-dependent Hamiltonian parameters as functions 
of rt = r{t) such as e{rt) and A(rt) (or as functions of the 
quantum state at the time). In this section, we present a 
few examples of the continuous measurement of a qubit 
with feedback, one with a linear feedback in the form 
A(rt) = Aq + Airt, and another with a state-dependent 
linear feedback to stabilize a qubit’s oscillation. 


A. Linear feedback, most likely path, and its phase 
space diagram 

Let us consider an instantaneous feedback loop in a 
qubit measurement introduced in section HI, with e = 0 
and the qubit Rabi frequency being a function of the 
measurement readout. We assume a linear form of the 
Rabi frequency as 

(52) 


where r* = r{t) is the readout as a function of time and 
Ao,Ai are constant parameters characterizing the bare 
Rabi oscillation and the linear feedback, respectively. We 
will see later in this subsection that this type of feedback 
can stabilize pre-determined (arbitrary) quantum states. 

This feedback loop in Eq. (52) has the advantage of 
fast processing because the qubit Hamiltonian depends 
directly on the value of the measurement readout, requir¬ 
ing no knowledge of the qubit state. For simplicity, we 
limit our discussion to an ideal case in which the qubit is 
measured with an efficient detector and no extra environ¬ 
ment dephasing. In this regime, we can re-parametrize 
the qubit Bloch vector into a single state parameter 9 
where z = cosd and y = sin0 (we set a; = 0 in this case). 
We note that, with this linear feedback Eq. (52), a typi¬ 
cal step size for the qubit state defined as A6 has average 
drift (Ad) ~ St and diffusion ((Ad)^) ~ St of the first or¬ 
der of the time step. Thus, the step size is still small. 

The statistics of quantum trajectories with this in¬ 
stantaneous feedback can be investigated following the 
same procedure presented in previous sections for the no¬ 
feedback qubit measurement. We start with writing the 
joint probability distribution of the trajectories, and then 
transform it into a path integral form. The action of the 
stochastic path integral with the feedback Rabi frequency 
^fb{rt) is given by, 

5/6 = ^ dt{-pe(d-Ao-Ai rH-r sind/r„) 

- (r^ - 2r cosd + l)/2Tm}, (53) 

where pe, d, and r are functions of time, omitting the 
time argument. This is a generalization of the action 
given in Ref. [18] where the quantum jump in qubit mea¬ 
surement is analyzed. We note here that in Eq. (53), 
we have chosen to use the (5t-expansion in both the state 
update and the readout probability distribution, because 
we are interested in investigating the action-extremized 
solutions, the optimal paths of the system. 

The optimal paths of the action in Eq. (53) are ob¬ 
tained by extremizing the action over all variables pg, d, 
and r, leading to two ordinary differential equations and 
one constraint, 

d = Aq + Ai r - r sind/rm, (54a) 

pe =percosd/rm + rsind/Tm, (54b) 

r =peAiT„,-pe sind + cosd, (54c) 

with arbitrary boundary conditions on the variable d and 
Pe. The solutions of these equations are the most likely 
paths of the system. By writing the action in the form 
Sfb = /dt(-ped +'H[pe,d,r]), where "H =pe(Ao +Air- 
r sin OlTra) - ("c^ - 2 r cosd + l)/2Tm is a stochastic Hamil¬ 
tonian, we can examine the most likely paths by using 
phase space analysis [18], motivated by the phase space 
concept in classical mechanics. Substituting the r con¬ 
straint in Eqs. (54c) into the action (equivalent to inte¬ 
grating the path integral over the variable r(t)), we ob¬ 
tain the action and the stochastic Hamiltonian in terms 


Afb{rt) = Ao + Airt, 
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d-Pe phase space Trajectories 0(t) and most-likely paths 



FIG. 4. (Color online) The phase space portrait and the most likely paths for the continuous measurement with linear feedback, 
(a) The 9-pe phase space portrait for the pure linear feedback case (Aq = 0) with Ai = 0.8 where the collapse angles shown 
as two vertical lines are at 6 si,s 2 ~ O.StTjO.Ttt. The curves represent the most likely paths for different stochastic energies: 
E = Q (long dashed blue), E = = -l/2rm = -0.5 (solid black), Ec < E < 0 (short dashed red), and E > 0 (dotted gray). 

The arrows show directions of the quantum state evolution along the paths. These most likely curves have the same form as 
in the plain measurement case, but with tunable collapse points, (b) Samples of numerically simulated qubit trajectories are 
plotted (solid gray fluctuating curves) along with the most likely paths for three different initial states: Oj = O.lrr (solid green), 
9i = O.Stt (dashed purple), and 9i = 1.27r (dotted orange), and the collapse states are shown as thin dashed black horizontal 
lines. These most likely paths are computed numerically from Eqs. (54), using the initial states 9{t = 0) = 0/ and the final 
condition pe(T) = 0 (the most likely paths without final fixed quantum states). For each set of boundary conditions (the 3 
sets), multiple solutions are found. Both of the dotted orange curves represent local maximum where the bottom one is more 
likely than the other. For the solid purple set, both curves are equally likely. For the dashed green curve, the initial condition is 
close enough to an attractor that we only find one most likely solution. All curves in the three sets correspond to the stochastic 
energies Ec < E < 0. 


of only the system variable 0 and its conjugate pg, 


{pgAiTjn+COs9-pgSm0)'^ 1 

n= --- +pgAo-—. (55) 

^ ' m. ^ > m. 


This quantity is explicitly time-independent and so it is 
a constant of motion for the most likely path, a solution 
of the ordinary differential equations in Eqs. (54). Let us 
define the stochastic energy E = 71, and then solve for the 
conjugate variable pg{0, E) as a function of 9 and E from 
Eq. (55). Each value of E will then correspond to a curve 
in the phase space portrait, describing dynamics of the 
most likely paths, the same way the individual classical 
trajectories are depicted as constant-energy curves on a 
phase space plot. 

As an example, we show in Figure 4(a) the phase space 
portrait for the pure linear feedback case (i.e., Aq = 0), 
where we plot the conjugate variable pg as a function of 

e, 


Pe{d,E) 


- cos 9 ± \/l + 2ETm 
AiTm - sin 0 


(56) 


for different values of E. This is an interesting case. 
We can see from the phase space plot in Figure 4(a) that 
there exists some attractors to which all states eventually 
limit to. 


These attractors coincide with the divergence of the 
conjugate variable pe{9,E) in Eq. (56), and also appear 
as stationary points where 0 = 0 in Eqs. (54a). Solving 
these equations, we obtain the attractors are located at 
(^si = j 7 r + arcsin(AiTm) and 9s2 = (j + l) 7 r-arcsin(AiTj„) 
where j = 0,2,4,... is an even integer. These attractors 
may be interpreted as stabilized states achieved by turn¬ 
ing on a linear feedback 0 < lAir^l < 1 and turning off 
the bare qubit frequency Aq = 0. They are effectively 
the new collapse points, rather than at the poles of the 
Bloch sphere, i.e., 9 = 0,7r,27r,... {z = ±1). Considering 
only angles between 0 and 27r, for the positive feedback 
0 < AiTm < 1, as jAiTml increases from zero, the stabi¬ 
lized states 9si, 0 s 2 move from 0 and tt toward each other 
and coalesce at 9si = 9s2 = 7r/2 when Air^ = 1; whereas 
in the negative feedback -1 < Ait^ < 0 , as |AiTm| grows, 
the stabilized states move toward each other in the re¬ 
gion of TT and 27r and coalesce at 3tt/2 when Air^ = - 1 . 
We simulate numerical qubit trajectories using the Monte 
Carlo method (Appendix E), showing that individual 
trajectories initialized at different states are eventually 
pinned to the stabilized states as predicted from the most 
likely path phase space. This is shown in Figure 4(b) for 
three different initial conditions, along with their most 
likely paths. 
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B. Stabilizing Rabi oscillation and its correlation 
function 


In this subsection, we show an example using the path 
integral to compute a correlation function for a system 
with linear feedback loop. The feedback loop of this ex¬ 
ample resembles the one used in the solid-state qubit 
measurement in Ref. [13] which has been adapted and 
realized in the transmon qubit experiment [14]. The the¬ 
oretical setup is an oscillating qubit (a double quantum 
dot) continuously monitored by a near quantum-limited 
detector (a quantum point contact), of which the Hamil¬ 
tonian and the measurement operator are in the same 
form as presented in Section III. We as before consider 
the ideal symmetric qubit case where e = 0, and the qubit 
is measured with an efficient detector with no extra en¬ 
vironment dephasing. The qubit state is represented by 
a single parameter 9 where = cosd and y = sin0. The 
stochastic master equation for the qubit state, in the St- 
expansion similar to Eqs. (10) (Stratonovich form, as dis¬ 
cussed in Ref. [13]), is given by. 


9 = Afb{9) - sm9cos9/Tm sin^/r,)/^, (57) 

where 9 and ^ are functions of time and A fh is the feed¬ 
back Rabi frequency, which is assumed to be a function 
of the qubit state 9. We later neglect the second term on 
the right side of Eq. (57) because we consider a diffusive 
Rabi limit A » . 

The feedback protocol considered in this subsection 
(and also in Ref [13]) is a linear feedback designed to sta¬ 
bilize the quantum oscillation of the qubit state, against 
the random phase kicks due to the measurement. The 
desired qubit evolution is described hy y = sin{Aj_t) and 
z = cos{Adt) where we define A^ as the target oscillation 
frequency. The difference between the actual phase 9 and 
the target phase A^t, denoted as 69{t) = 9{t) - Adt, is 
used to control the oscillating part of the qubit Hamil¬ 
tonian. Therefore, we write the feedback Rabi frequency 
as, Afb = Ad{l-F59), where F is the dimensionless feed¬ 
back factor. This feedback loop continuously corrects the 
random changes of the state made by the measurement 
so that the outcome trajectory closely follows the desired 
qubit oscillation. 

Let us assume that the phase difference (5d is a slowly 
changing variable as compared to the oscillation with the 
desired frequency A^. Therefore, we can average the 
fluctuating process described in Eq. (57) over the oscilla¬ 
tion period 5t = 2 t:IA d, taking 59 to be constant during 
this period. The period 5t will eventually be our new 
time scale. We define a new noise ^ as a time-average 
of the last term of Eq. (57) over the oscillation period, 

|(t) = -{list) sm{59{t') + Adt')^{t')lTll^ with 

its zero ensemble average {^{t)) = 0 and its variance 
being (f(t)^) = (I/5t)^ sm'^{59{t') + Adt')lTm = 
ll2Tm5t. We then can simplify the differential equation 


Eq. (57) to one in terms of the phase difference 59, 


59{t)=-FAd59{t)+i{t), (58) 


with the new time scale 5t. 

Here we will use the above differential equation 
Eq. (58) to compute a correlation function Kz{t) = 
{z{t)z{t + t)) using the path integral approach. Follow¬ 
ing the derivation of the action in Eq. (3) and (8) where 
q = 59 and P(f) = (rmStln)^^^ exp{-Tmi^5t) are now 
our new system variable and a noise probability density 
function, we obtain the action in this form, 

Sfb = dt'{ - ips9{59 + FAo59 -1) - (59) 

omitting the boundary terms. Note that we have written 
the pure imaginary conjugate variable {ipse) explicitly 
with i. As before, we can compute an average quantity 
by integrating the paths (•••) = J\fjV59Vpsg'D^e^f'’---. So, 
we first write the correlation function Kz{t) in terms of 
the phase difference 59, then average over the oscillation 
period 5t, getting rid of the fast fluctuating parts, and we 
are left with ^^(r) « {cos[59{t)-59{t + T)])cos{AdT)/2 + 
sin[(I0(t) - 59{t + t)]) sm{AdT)/2. The correlation func¬ 
tion apparently can be written in terms of the real part 
and imaginary part of the following quantity, 

(60) 


where we have used the same notations, such as f'D59 
for the integral over all possible paths 59{t), and Af for 
a normalized factor. 

The Gaussian integral over ^ in Eq. (60) is quite 
straightforward resulting in a bilinear term in pse, which 
then leads to another Gaussian integral of pse- As one 
would expect, these integrals generate another prefactor 
that cancels the normalized factor Af. Gonsequently, the 
last integral over 59 is left as = jV59 

where the exponent S' is given by, 

5'= rAt'{-T^{59 + FAd59)'^ 

JO 

+ i59 5{t'-1) - i59 5{t'-1 - t)}. (61) 

This effective action can be transformed further using the 
integration by parts, for example, fdt5959 = - fdt5959, 
assuming that the boundary conditions for 59 at both end 
points vanish. We then write the action in terms of the 
inverse Green’s function and the source term as discussed 
in section V A, S' = -I Jdt'dt"59{t')Ggl{t',t")59{t") + 
fdt'Jse{t')59{t') where the inverse Green’s function, the 
Green’s function, and a particular form of the source term 
are given by, 

G-sl{t',t") = 2Tm5{t' -1") {-dl, + F^Al), (62) 

Gse{t',t") = l/{dT^\FAd\)exp{-\FAd{t' -t")\} , (63) 
Jse{t') = i5{t' -t) - i5{t' -t-r). (64) 
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Performing the last Gaussian functional integral over 
59 gives, 




) =exp|- y~ dt'dt"J,t") J{t") 

e-|^A,|r _ ^ 


= exp 


4Tm\FAd\ 


(65) 


for the positive value of the time difference t. This quan¬ 
tity is a real number, therefore we obtain the correlation 
function, 


Kz{t) 


cos(AdT) 

-^-exp 


g-FAaT_ 

dTmFAd r 


( 66 ) 


assuming that FA^ > 0. This result agrees with the so¬ 
lution found in Ref. [11] which is presented in different 
notation. 


and we have computed the correlation function for the 
qubit trajectory using the path integration method. 

So far, the stochastic path integral formalism in the 
context of continuous quantum measurement has been 
proven to be useful in studying the statistics of quantum 
trajectories; however, there are some unsolved issues that 
need to be further explored. One is the limitation of the 
statistical average solutions derived from the perturba¬ 
tion expansion theory. Only the Hrst few orders of the 
expansion have been computed, resulting in the solutions 
that are valid only in the certain parameter regimes. We 
hope to Hnd solutions in an arbitrary regime, possibly 
with some modifications of our approach. Another issue 
is the assumption of the instantaneous feedback, which 
can be difficult to realize in experiments. Feedback loops 
modelled with time delays will be taken into account in 
future work. 


VII. CONCLUSION 

We have developed and extended the stochastic path 
integral technique to study statistical behaviour of a 
quantum system under weak continuous measurement, 
as well as measurement with feedback, presented with 
several qubit examples. The path integral approach is 
constructed based on the joint probability distribution 
of the measurement records, which is then extended to 
the distribution of quantum states, describing all possi¬ 
ble quantum trajectories. We have shown that with this 
path integral and its action formalism, the optimal dy¬ 
namics, such as the most likely paths, can be obtained 
naturally from the extremization of the action, whereas 
other statistical quantities can be achieved from direct 
integration or perturbation theory. In the case of plain 
measurement of a qubit, we have derived analytic so¬ 
lutions for the average trajectory, the variance, and the 
correlation functions conditioning on the Hxed initial and 
final states, which show an excellent agreement with the 
numerically simulated data. 

We have also presented a diagrammatic perturbation 
method used in computing expectation values and cor¬ 
relation functions of quantum trajectories, and elabo¬ 
rated it with examples of the qubit with Rabi oscilla¬ 
tion case. The variances and multi-time correlation func¬ 
tions of qubit trajectories in the short-time regime have 
been revealed using this method given initial conditions, 
and the results are in good agreement with the numer¬ 
ical simulation. Moreover, we have considered quan¬ 
tum measurement with feedback control, using the ac¬ 
tion principle to investigate the dynamics of the most 
likely paths of a qubit with linear feedback on its os¬ 
cillating frequency. We have discovered that the direct 
linear feedback, manipulating the qubit Hamiltonian in¬ 
stantly using the measurement readout, can stabilize the 
qubit state to arbitrarily chosen pure states. We have 
also considered the example of the feedback loop stabi¬ 
lizing the qubit Rabi frequency introduced in Ref. [13], 
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Appendix A: Connection to other path integral 
formalisms for quantum measurement 

There are interesting comparisons between the path 
integral formalism we introduced in the main text and 
other path integral approaches built upon Feynman path 
integral in quantum mechanics [46]. As we have shown, 
our formalism is aimed at describing the probability dis¬ 
tribution of quantum trajectories, paths of state of a sys¬ 
tem under continuous measurement, on its quantum state 
space (such as Hilbert space for pure states). Each in¬ 
dividual quantum trajectory in the path integral is re¬ 
alizable and tractable, as demonstrated experimentally, 
such as, in solid-state systems [17, 19]. However, for 
other path integral formalisms developed from the Feyn¬ 
man path integral to investigate quantum systems under 
measurement [22-24, 26, 47], the evolution of the quan¬ 
tum state (or wavefunction) is based on interference of all 
possible classical paths in measurable configuration space 
(such as the system’s positions or spin states). Thus, in 
this latter case, the integrations are over conhguration 
coordinates of the measured system. 

Despite the differences in forms, the two approaches 
mentioned above can be related. The path integral via 
Feynman’s concept can be used to compute the proba¬ 
bility distribution of the measurement results, which is 
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then, as shown in the main text, one of the most im¬ 
portant ingredients in constructing our stochastic paths 
in quantum state space. To elaborate this connection in 
more detail, we consider the path integral method pre¬ 
sented by Caves [24], for measurements providing infor¬ 
mation about the position x{t) of a nonrelativistic, one¬ 
dimensional quantum system, evolved in time. In that 
approach, the effect of measurements is to restrict the 
sum over paths by weighting each path differently de¬ 
pending on the measurement results. These weights ap¬ 
pear in the path integral in Ref. [24] as the ‘resolution 
amplitude’, T{x{t) - x{t)), which also accounts for the 
imprecision of the measurements. 

Let us assume instantaneous position measurements 
that are equally distributed in times, at to,ti, 
where tk = to + kSt, giving measurement readouts denoted 
by xo,xi, ...,Xn-i- The joint probability amplitude $ of 
the measurement records, and that the system is at Xn at 
time tm given an initial wavefunction is in this form. 




ilV’o) = T(a;fc -Xfe)j 


n-1 , \ 

n (a;fc+i[e"*'^‘-^lxfe) (xqIV’o), (A1) 

fc=o / 


where the position coordinates are denoted by 
xo,xi, ...,Xn at time and the integral 

measure is defined as J'Dx{t) = f da;o-"da:„_i. We 
note here that we have modified the ordering of the 
measurements from the original version by Caves, in 
which he assumed that the measurements start only 
after the initial state has evolved for time 6t. We 
assume that this change has an infinitesimal effect on 
the probability amplitude as St ^ 0. The first bracketed 
term on the right hand side of Eq. (Al) describes the 
influence of the measurement. Without it, the path 
integral will be the usual Feynman path integral, exactly 
equal to {xn\e~^’^^\'ipo) or the wave function at the final 
time tn = T. 

From Eq. (Al), the joint probability distribution func¬ 
tion (PDF) of the sequence of measurement readouts can 
be obtained by integrating the square of the probability 
amplitude over the final coordinate 

P{xo,Xi,...,Xn-l\'P’o) = dXn\^{xo, ■■■,Xn-l]Xn,tn\lpo)f , 

(A2) 


given the initial wavefunction 'tpQ. 

In order to see that this joint PDF in Eq. (A2) is the 
same as what we have earlier in the main text (in this spe¬ 
cial case), we first rewrite this term /da:o(a;i|e"*'^‘'^T(a;o- 
a^o)|a^o)(3^olV’o) = (xilUstMoltpo) as a probability ampli¬ 
tude right after one measurement (Mq) and one uni¬ 
tary transformation (Ust)- Then, we obtain a probabil¬ 
ity density function for the first measurement outcome 
as P(a;o|V’o) = /dxiKxilD^tMolV’o)!^- We further intro¬ 
duce an updated wavefunction \pk+i), a result of a state 


[ipk) that has gone through one measurement and one 
unitary transformation Usti 


IV’fc+i) 


/ <lxke~'^^*^T{xk - xk)\xk){xk\ipk) 
\/P(Xfe]V’fc) 

U5tMk\'ipk) 

\Jp{xk\i’k) 


where the denominator assures the correct norm of 
the wavefunction. Then, the joint distribution of 
the measurement readouts of the first two times is 
given by P(xi,xo[f/'o) = jdx 2 \{x 2 \UstMi\'p-i)\^P{xQ\'po) = 
P(a;i]'i/'i)P(a;o|V’o)- Repeating this procedure further to 
the next measurement readouts X 2 ,X 3 , ...Xn-i, at the end, 
we obtain the full joint PDF, 


n-1 

P(xo,Xi, ...,Xn-lllpo) = Yl 


k=0 


(A4) 


as a product of the conditional probability density 
functions, defined in a general form as Pixkltpk) = 
fdxk+i\{xk+i\UstMk\p’k)f- This is the joint PDF of the 
measurement readouts, the main component of the joint 
PDF in Eq. (1). 

To obtain an analogous form of the joint PDF in 
Eq. (1), we add probability density functions of the quan¬ 
tum state (or wavefunction) as delta functions to every 
single time step. We obtain the full joint PDF, 


n-1 

'P^p= Yl P{i’k+l\'lpk,Xk)Pixk\lpk), (A5) 

k=0 

where the update state in this case is written as 
P{'ipk+i\'ipk,Xk) = S^ll-ipk+i) - UstMk\'ipk)/\/P{xk\tpk)} as 
in Eq. (A3). The integer d is the dimension of the vec¬ 
tor \^pk) (which maybe generalized to infinite dimension 
via a functional form of the 5-function), and is called 
a joint probability density function of the measurement 
outcomes {xk} and the wavefunctions {jV’fc)} of the mea¬ 
sured system. 

We note that our path integral presented in the main 
text deals with mixed states instead of pure states, and 
the system we consider is the qubit (or spin) system 
with the discrete basis, i.e., / dxla;)(a;l ^ Es=i,-i |s)(s|. 
Therefore, the operator Mk in the discussion above is 
equivalent to what we have as the measurement operator 
Mk = Mrk = T(rfc - I)|l)(l| + T(rfc + I)] - I)(-I|, defined 
in Section III, and the resolution function is a Gaussian 
function, T(rfe-s) = (5t/27rTm)^'^‘^exp{-5t(rfe-s)^/4rm} 
for s = 1 and s = -1. 

It is also worth mentioning that in our approach, the 
stochastic trajectories in quantum state spaces can be 
thought of as classical trajectories in configuration space, 
such as trajectories on the Bloch sphere can be con¬ 
sidered as three-dimensional random walks in a unit- 
radius sphere. Then, one could find connections between 
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our path integral and the formalism in classical stochas¬ 
tic processes, such as the Pilgram-Sukhorukov-Jordan- 
Biittiker path integral [29], the Martin-Siggia-Rose for¬ 
malism, the Wiener integral, or Feynman-Kac path inte¬ 
gral (see Eqs. (16),(17), and (18)), however, we do not 
cover the discussion in this paper. 


Appendix B: The path integral’s action for the qubit 
measurement 

We show the derivation here, how the action in Eq. (17) 
can be written in this form, 

S[u\ = S[u\-^ (Bl) 

where u is the optimal path, extremizing the action iS[u]. 

To see this, one can Taylor expand the action Eq. (17) 
in discrete forms, such as 5 [m] = S[u+r]] = 5['u]+5'['«]?7 + 
^S'[u]ri'^ + and show that higher order terms van¬ 

ish {0{rf) = 0). However, there is a simpler way to see 
this, using the action in a continuous form, 

S = - f -ta.nhu(t)u(t) +—^\, (B2) 

Jo [ 2 2Tm) 

which is the same action as in Eq. (17), with definitions of 
the time integral, 5tYjk=o = Jq and the deriva¬ 

tive, {Ak+i-Ak)ISt = A{t). In this form, one can see that 
the last two terms in the action can be integrated easily. 
The second term can be written in this general form, 

f dt = f dtF{u{t)), 

Jo Jo 

= F{u{T))-F{u{0)), (B3) 

where F{u{t)) = /dutanhu = ln(coshM(t)). The 

contribution of this term to the action is only de¬ 
pendent on the boundary terms. Therefore, we can 
write the last two terms in Eq. (B2) as equivalent to 
f^dt {ta.TLhu{t)u{t) - l/2Tjn}- 

After this simplihcation, the first term in Eq. (B2) can 
be written as - /q dt^ + 2.u{t)'f]{t) + ?7(t)^}, per¬ 

forming the integration by parts giving J^dtu{t)fi{t) = 
- j^dtu{t)ri{t) = 0. The action is then given by, 

S[u\ = S[u\-'^jdt'q{tf, (B4) 

where, 

iS[u] = - f dt I- ta.nhu{t)u{t) + -—|, (B5) 
Jo ( 2 2TmJ 

is the action in terms of the optimal solution u{t). 


Appendix C: Inverse matrix M ^ 


We show here the inverse of the matrix M in Eq. (21), 
assuming that it is a square matrix with d dimensions. 
The inverse matrix is given by, 


( d d-1 

d-2 - 

1 

d-1 2{d-l) 

2{d-2) - 

2(d-d+l) 

d-2 2{d-2) 

3(d-2) - 

3 

, 1 2 

3 

d , 


where the diagonal elements are of the form = 

{dtlTm)k{d - {k - 1))/{d + 1) = {5tjTm)k{n - k)ln, where 
d = n - 1, and the other off-diagonal elements are M~^ = 
^kj = - k)/n for k > j as verihed by direct 

calculation. 


Appendix D: Full solutions for the qubit 
measurement without Rabi oscillation 

The solutions of the preselected and postselected av¬ 
erage and variance in z to infinite order are shown here. 
The average is given in terms of the infinite summation 
as. 


Zp 




E 


m=0 


1 

(2m)! 


j2m 


duf^ 

\ J 


tanh(uj) 



where the differentials in the bracket are evaluated at the 
optimal path u. Similarly, the variance is given by, 

Zp {Zj ) Zj~ Zp{zj) Zj 


= E 


1 


j2m 


(2m)! du‘j'^ 

-It 1 / 

\r^0 (2w)! 


tanh^{uj) 


tanh(Mj) 


(-if)) ■ 

(D2) 


The average are given explicitly in Eq. (27) in the 

main text. 


Appendix E: Numerical simulation 

The numerical data presented in Figure 1, 3 and 4 
are from the simulation of quantum trajectories using 
Monte Carlo method. The numerical trajectories are gen¬ 
erated in n-discrete steps of St. Starting from an initial 
state go7 each step, we randomly generate a measure¬ 
ment readout from a distribution, for example, P{rk\qk), 
and compute a quantum state from the update equation 
Qk+i = S[qk,rk] (or pk+\ = O^UstMrdPk]), repeating 
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the procedures from /c = 0 to n - 1 to get a full trajectory 
{g/clo- 

In Figure 1, we simulate the data using n = 100 time 
steps, with the step size 6t = 0.006 and = 1- In Fig¬ 
ure 3, we use n = 500, St = 0.01, and = 10, where the 
Rabi frequency A is set to 27r. We note that in the main 
text we present these numbers in the unit of Tm- 

For the data presented in Figure 4, where the Rabi os¬ 
cillation is linearly dependent on the highly fluctuating 
measurement readouts, the state update equation needs 
to be modified to minimize numerical errors in each time 
step of the calculation. The update state is then com¬ 
puted from pk+i where 

the measurement operator and unitary operator Ust 
are each divided into 10 pieces and are operated onto the 
qubit state pk alternately. We simulate the data using 
n = 1000 time steps, with the step size St = 0.01 and 
Tm = 1. The feedback Rabi oscillation is A{tk) = O.Sr^, 
where Vk is a generated measurement readout at time tk- 


Appendix F: White noise limit 


that the distribution of rjk is Gaussian with the same 
variance TmISt and the mean Zk- 

Appendix G: The initial-boundary source terms 

Let us assume that x is the only system variable we 
consider, where an update equation is written as Xi+i = 
Xi -H L{xi,^i)St and the probability density function for 
a measurement readout is proportional to We 

then write the action for this system as, 

n-l r 1 

^ " S \ -Pi{xi+i -^PiL{xi,^i)St- -^iSti. (Gl) 
i=0 ^ ^ J 

where p is the conjugate variable. In this case, 
we could write the first term, Y,7=o ~ ^i) - 

-Er=o^ 'L]=oPzCijXj, where = Si+ij - Spj, however, 
the matrix C would not be a square matrix, and its in¬ 
verse is not simply defined. Therefore, we separate out 
one term, pqXq, from the sum, which symmetrizes the 
double sum. 


In the main text, the state updating procedure is based 
on the assumption that a measurement outcome is ran¬ 
domly generated from a readout distribution P{rk\qk), 
which is a function of a state right before the measure¬ 
ment. However, in the limit when the measurement read¬ 
out is highly fluctuating, for example, the standard devi¬ 
ation of the readout distribution is much larger than the 
measurement response, (TmlSt)^^"^ » 1, we can approx¬ 
imate the readout as a sum of two parts, one being its 
average (which is related to the measured system state) 
and another being a zero-mean independent fluctuating 
noise. Following the readout distribution in Eq. (9), the 
average of the readout at a time tk is exactly (r^) = Zk, 
therefore, we write Vk = Zk + rjk where rjk is a zero-mean 
Gaussian white noise. 

The next step is to find a probability distribution for 
the independent noise rjk, which we then approximate 
as having the same variance as the original distribution, 
which in this case the variance is TmjSt. This approx¬ 
imation is valid as long as the variance is much larger 
than the separation between two outcomes r = 1 and 
r = -1. A more rigorous proof can be done by writ¬ 
ing Pk = \/^^k where ^ is a Gaussian white noise with 
variance St~"^. Note that ^ is the time-derivative of the 
Wiener increment, scaling as SP^!"^ (ltd calculus) [42], 
therefore, in any expansion, we need to keep terms con¬ 
taining r^St^ Ri TmSt. Using these rules, we can write, 

P{rk\zk)=MiP^e-^^^^-^^" + 

RiA/'exp|-^(rfc-2fc)^|, (FI) 

where these two equations are exactly equal in the expan¬ 
sion up to the first order of St. The second line confirms 


n-l n-l n 

Y, -Pi{Xi+l - Xi) = PoXo - Y (G2) 

2=0 2=0 j = l 

resulting in a square matrix (G~^), 


(G -0 




/I 0 0 
-1 1 0 
0 -1 1 


0 


\ 0 - 0 -1 1 / 


(G3) 


with the row index, f = 0 to n -1, and the column index, 
j = 1 to n. For the measurement readout term, we define 
its square matrix. 


^St 


(1 

0 

0 


0 0 
1 0 
0 1 


0 


\0 - 0 0 1 / 


(G4) 


with the row and column index being i, j = 0 to n - 1. 

As a result, we can rewrite the action S into two sep¬ 
arated terms, a free action Sp and an interaction action 
Si, where 


n-l n 1 n-ln-l 

= - E Zp^(G-\^x, --YZ (G5) 

i=0 j = l ^ i=0 j=0 

PoXoSip + piL{x^,^i)St 

i=0 t 

The appearance of the first term in the interaction ac¬ 
tion, pQXgSifi, is the reason we have the extra terms (we 
represented them as B) in Eq. (41b). 
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Appendix H: The left continuous Heaviside step 
function 


We now show why the two-point correlation functions 
derived from our path integrals, such as the one derived 
from the free action in Eq. (G5), that is 


{x{t)p{t'))F = (HI) 


consist of a left continuous Heaviside step functions 0(t) 
that behaves differently from the usual Heaviside step 
function. This left continuous Heaviside step function 
has the properties, 0 ( 0 ) = 0 and limi^o+ 0 (^) = 1 ) 
while the usual Heaviside step function has properties, 
0(0) = 1/2, limt^o+ 0(0 = 1) and limj^o- 0(0 = 0- This 
is true from a point of view of the discrete form of a cor¬ 
relation function, {xaPb)F = {Gx)a,b, that the propagator 
has this property {Gx)a,b = 1 when a> b and it vanishes 
otherwise. 

Let us start from writing a free generating function in 
a discrete form, Zf[Jx, Jp, J.c] 


/ 5t ^ ^ ^ ^ f Tl _1 7T, 


i=0 j=l 


n-1 n I 

fe=0 k=l 


^ [ -| n-1 n-1 n-1 

/ d[Cfe]o"^expi--^ Z O(G^0i.iO + Z • 

■J I ^ i=0 j=0 


(H2) 


k=0 


These matrix integrals can be carried out using multi¬ 
dimensional Gaussian integrals. The first one is. 


J' dpdx 


^-p Ax+p Jp + J^x _ A Vf 


(27r0^ 
Det[H] ’ 


(H3) 


where p is a pure imaginary d-dimensional vector, 
x, Jp, Jx, are real vectors, and A is a real and invertible 
matrix. Another is. 


fd(e-i 


A = (H4) 


(27r)^/2 


(Det[H])i/2’ 


where ^ is a real d-dimensional vector and H is a real and 
symmetric matrix. After integrating over all variables, 
the free generating function is given by, 

i n n-1 I 

E 'EiJx)iiGx)^,jiJp)jSf [ 

i=ii=o j 

( 1 n-1 n-1 I 

, E E(T?)*(Qk.(Tc),dt" . (H5) 

^ i=0 j=o j 

We note that the integrals generate a prefactor 
Det^G^i) ^ (^TTi)", knowing that Det(G”^) = 1, and an¬ 


other prefactor ) ^, which both are then canceled with 
the prefactor in Eq. (H2). Also, from Eqs. (G3)-(G4), we 
can compute the inverses of them as, 


{Gx)i,j - 


/I 0 0 - 0\ 
110-0 
111-0 

\i - 1 1 IJ 


(H 6 ) 


with the row index, f = 1 to n, and the column index, 
j = 0 to n - 1 , and 


(Q)*o = Yt 


/I 0 0 - 0\ 
0 10-0 
0 0 1-0 

(0 - 0 0 1/ 


(H7) 


with the row and column indices, i, j = 0 to n - 1 . 

Now, let us compute the two-point correlation function 
{XaPb)F, 
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^ [ n-1 n 1 n-1 n-1 

(a;aP6)F sAA/d[a:fc]"d[pfc]o"M[Cfc]rHa^aP&)exp]- ^ '^Pi{Gl^)i,jXj - - Y, E 

i=0 1 = 1 ^ i=0 1=0 


1 a 1 a 


Z F \^Jx^ Jp^ 'A^] 


d{Jx)a St d{Jp)b 

1 a 1 a 

St d{Jx)a St d{Jp)b 

[1 , if, a > & + 1, or, a > 6, 


exp 


Jx tJp “0 

n n-1 

E '^i’^x)i{Gx)i,j{Jp)jSt 

Li=l j=0 


-{Gx)a,b - 


0 , otherwise. 


Jx,Jp=0 


(H8) 


where Af = (|^)^ /(27i0”- The conclusion in the last 
line comes from the fact that the square matrix, (Gx), 
has the row and column indices different by one time step 
(see Eq. (H6)), which actually originates from the indices 
of, {xk}i, and {pk}o~^- 


Appendix I: The 10 diagrams of the correlation 
function (w(ti)w(t 2 )) 

Continued from Eq. (51), we presented all 10 diagrams 
here. 


{v{ti)v{t2))^°^ = (^)(^l)^)(^2)e‘^0F^ 



-vf Gy{tl,0)Gy{t2, y' ^Gy{tl,t')Gy{t2,t') { ^2 + UV^J G y{t' , O)^ + aVjWjGy{t' , 0)Gyj{t' , O)} 




fdt’ 




1^2-^ I 






+ aviwje ^ e 


( 11 ) 


where the definitions of the parameters vi,a,K 2 , ■■■ are 
discussed in Eqs. (38)-(39). 
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